gc()
## used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 409742 21.9 841634 45 NA 658114 35.2
## Vcells 793965 6.1 8388608 64 102400 1802080 13.8
rm(list=ls())
library(gam)
## Loading required package: splines
## Loading required package: foreach
## Loaded gam 1.20
library(scran)
## Loading required package: SingleCellExperiment
## Loading required package: SummarizedExperiment
## Loading required package: GenomicRanges
## Loading required package: stats4
## Loading required package: BiocGenerics
## Loading required package: parallel
##
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, basename, cbind, colnames,
## dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
## grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
## order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
## rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
## union, unique, unsplit, which, which.max, which.min
## Loading required package: S4Vectors
##
## Attaching package: 'S4Vectors'
## The following object is masked from 'package:base':
##
## expand.grid
## Loading required package: IRanges
## Loading required package: GenomeInfoDb
## Loading required package: Biobase
## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.
## Loading required package: DelayedArray
## Loading required package: matrixStats
##
## Attaching package: 'matrixStats'
## The following objects are masked from 'package:Biobase':
##
## anyMissing, rowMedians
##
## Attaching package: 'DelayedArray'
## The following objects are masked from 'package:matrixStats':
##
## colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges
## The following objects are masked from 'package:base':
##
## aperm, apply, rowsum
library(ggplot2)
library(dynamicTreeCut)
library(umap)
library(ComplexHeatmap)
## Loading required package: grid
## ========================================
## ComplexHeatmap version 2.4.3
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
##
## If you use it in published research, please cite:
## Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional
## genomic data. Bioinformatics 2016.
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================
library(gridExtra)
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:Biobase':
##
## combine
## The following object is masked from 'package:BiocGenerics':
##
## combine
library(grid)
library(viridis)
## Loading required package: viridisLite
library(destiny)
##
## Attaching package: 'destiny'
## The following object is masked from 'package:SummarizedExperiment':
##
## distance
## The following object is masked from 'package:GenomicRanges':
##
## distance
## The following object is masked from 'package:IRanges':
##
## distance
library(circlize)
## ========================================
## circlize version 0.4.12
## CRAN page: https://cran.r-project.org/package=circlize
## Github page: https://github.com/jokergoo/circlize
## Documentation: https://jokergoo.github.io/circlize_book/book/
##
## If you use it in published research, please cite:
## Gu, Z. circlize implements and enhances circular visualization
## in R. Bioinformatics 2014.
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(circlize))
## ========================================
library(ggforce)
library(edgeR)
## Loading required package: limma
##
## Attaching package: 'limma'
## The following object is masked from 'package:BiocGenerics':
##
## plotMA
##
## Attaching package: 'edgeR'
## The following object is masked from 'package:SingleCellExperiment':
##
## cpm
library(fpc)
library(Seurat)
## Attaching SeuratObject
##
## Attaching package: 'Seurat'
## The following object is masked from 'package:SummarizedExperiment':
##
## Assays
library(scater)
##
## Attaching package: 'scater'
## The following object is masked from 'package:limma':
##
## plotMDS
#base='/Users/gabriele.lubatti/Desktop/Phd/Cell_Competition'
base_path="/Users/gabriele.lubatti/Desktop/Phd/Cell_Competition/Cell_Competition_Github/Input_data_Paper_Figures"
#setwd("/Users/gabriele.lubatti/Desktop/Phd/Cell_Competition/Script")
#load("data_tpm.Rda")
#load(file="data.Rda")
#load(file="good_cells.Rda")
#load(file="meta_dati.Rda")
#geni_comuni=intersect(row.names(data),row.names(data_tpm))
#data_raw=data[geni_comuni,colnames(data_tpm)]
#data_tpm=data_tpm[geni_comuni,]
setwd(base_path)
load(file='data_raw.Rda')
load(file='data_tpm.Rda')
load(file='QC_start.Rda')
load(file='meta_dati.Rda')
load("HVGC_giusto_no_un.Rda")
load("Gene_Name_ID.Rda")
data_p53=read.csv("p53_full.csv",sep=";",header=T)
up_p53_mouse=read.csv("up_p53_topo.csv",sep=",",header=T)
down_p53_mouse=read.csv("down_p53_topo.csv",sep=",",header=T)
geni_upr=read.csv2("ana_upr.csv",header=F)
#setwd("/Users/gabriele.lubatti/Desktop/Phd/Cell_Competition/SC_course")
#load(file="data_tpm_course.Rda")
#load(file="data_raw_course.Rda")
#load("meta_dati1.Rda")
#load("HVGC_giusto_no_un.Rda")
#load("Gene_Name_ID.Rda")
meta_dati=meta_dati[colnames(data_raw),]
Condition=rep(0,length(meta_dati$Sample))
Condition[grep("FMK",meta_dati$Sample)]="Cell competition OFF"
Condition[grep("DMSO",meta_dati$Sample)]="Cell competition ON"
Extended Figure 1 A
QC_start_good=QC_start
QC_start_good$quality_control="Good Cells"
pw1=ggplot(QC_start_good, aes(QC_start_good$quality_control,QC_start_good$numAbove10rpm)) +
geom_boxplot() +ggtitle("QC (723 cells)")+ylab("numAbove10rpm")+xlab(NULL)
pw2=ggplot(QC_start_good, aes(QC_start_good$quality_control,QC_start_good$fracSpikes)) +
geom_boxplot() +ggtitle("QC (723 cells)")+ylab("Fraction of spikes")+xlab(NULL)
pw3=ggplot(QC_start_good, aes(QC_start_good$quality_control,QC_start_good$fracMTs)) +
geom_boxplot() +ggtitle("QC (723 cells)")+ylab("Fraction of MTs")+xlab(NULL)
pw4=ggplot(QC_start_good, aes(QC_start_good$quality_control,QC_start_good$fracGenes)) +
geom_boxplot() +ggtitle("QC (723 cells)")+ylab("Fraction of Genes")+xlab(NULL)
pw5=ggplot(QC_start_good, aes(QC_start_good$quality_control,QC_start_good$fracMappedReads)) +
geom_boxplot() +ggtitle("QC (723 cells)")+ylab("Frac of Mapped Reads")+xlab(NULL)
pw6=ggplot(QC_start_good, aes(QC_start_good$quality_control,QC_start_good$totReadsLog10)) +
geom_boxplot() +ggtitle("QC (723 cells)")+ylab("Tot number of reads(Log10)")+xlab(NULL)
grid.arrange(pw1,pw2,pw3,pw4,pw5,pw6,ncol=3,nrow=2)
## Warning: Use of `QC_start_good$quality_control` is discouraged. Use
## `quality_control` instead.
## Warning: Use of `QC_start_good$numAbove10rpm` is discouraged. Use
## `numAbove10rpm` instead.
## Warning: Use of `QC_start_good$quality_control` is discouraged. Use
## `quality_control` instead.
## Warning: Use of `QC_start_good$fracSpikes` is discouraged. Use `fracSpikes`
## instead.
## Warning: Use of `QC_start_good$quality_control` is discouraged. Use
## `quality_control` instead.
## Warning: Use of `QC_start_good$fracMTs` is discouraged. Use `fracMTs` instead.
## Warning: Use of `QC_start_good$quality_control` is discouraged. Use
## `quality_control` instead.
## Warning: Use of `QC_start_good$fracGenes` is discouraged. Use `fracGenes`
## instead.
## Warning: Use of `QC_start_good$quality_control` is discouraged. Use
## `quality_control` instead.
## Warning: Use of `QC_start_good$fracMappedReads` is discouraged. Use
## `fracMappedReads` instead.
## Warning: Use of `QC_start_good$quality_control` is discouraged. Use
## `quality_control` instead.
## Warning: Use of `QC_start_good$totReadsLog10` is discouraged. Use
## `totReadsLog10` instead.

Figure 1 C
my.tree <- hclust(my.dist
, method="ward.D2" # met = "ward.D2" was adopted form the reference
)
my.clusters <- cutreeHybrid(my.tree
,distM=as.matrix(my.dist),deepSplit=0
,minClusterSize = 35)$label#5)$label
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
length(unique(my.clusters))
## [1] 5
classe<-data.frame(Cell_type = meta_dati$Sample
,cluster = my.clusters)
classe$condizione=Condition
table(classe$cluster,classe$condizione)
##
## Cell competition OFF Cell competition ON
## 1 78 229
## 2 77 94
## 3 125 9
## 4 71 0
## 5 16 24
Cluster=as.factor(classe$cluster)
sp<-ggplot(as.data.frame(object$layout),aes(umap_met[,1],umap_met[,2],color=Cluster)) + geom_point()
sp + scale_color_manual(values=c("#DD6400", "#800080","#0000FF", "#006400","#56B4E9","#A9A9A9"),labels=c("Normal (winner) Epiblast (78FMK/229DMSO) ","Extra Embrionic Ectoderm (77FMK/94DMSO)","Intermediate (125FMK/9DMSO)","Loser Epiblast (71FMK/0 DMSO)","Visceral Endoderm (16FMK/24DMSO)"))+labs(x = "Umap1",y="Umap2") + ggtitle(NULL)

#setwd("/Users/gabriele.lubatti/Desktop/Phd/entropy_of_mixing_library/teichman_data")
#load("gene_id_name_table.RData")
#row.names(temp)=as.character(temp$Gene_stable_ID)
#gene.name.id<-as.character(temp[,"Gene_stable_ID"])
#names(gene.name.id)<-as.character(temp[,"Gene_name"])
#geni_mt_id=row.names(data_raw)[row.names(data_raw)%in%as.vector(gene.name.id[grep("mt-",names(gene.na#me.id))])]
#cells_epi=colnames(data_raw)[classe$cluster==1|classe$cluster==3|classe$cluster==4]
#row.names(classe)=colnames(data_raw)
#classe_epi=classe[cells_epi,]
#classe_epi=classe_epi[classe_epi$condizione=="Cell competition OFF",]
#cluster_condition=rep(0,length(classe_epi$cluster))
#cluster_condition[classe_epi$cluster==1&classe_epi$condizione=="Cell competition OFF"]="Winner cells"
#cluster_condition[classe_epi$cluster==3&classe_epi$condizione=="Cell competition OFF"]="Intermediate #cells"
#cluster_condition[classe_epi$cluster==4&classe_epi$condizione=="Cell competition OFF"]="Loser cells"
#cluster_condition[classe_epi$cluster==1&classe_epi$condizione=="Cell competition ON"]="Win Epi DMSO"
#cluster_condition[classe_epi$cluster==3&classe_epi$condizione=="Cell competition ON"]="Int Epi DMSO"
#sum_mt=apply(data_raw[geni_mt_id,row.names(classe_epi)],2,sum)
#cluster_condition=factor(cluster_condition,levels=c("Winner cells","Intermediate cells","Loser #cells"))
#data_plot=data.frame(cluster_condition,sum_mt)
#p<-ggplot(data_plot, aes(x=cluster_condition, y=sum_mt, fill=cluster_condition)) +
#geom_boxplot()+labs(fill="Cluster")+ylab("Number of reads mapping to mtDNA")+xlab("Cluster ")
#p+scale_fill_manual(values=c("#DD6400", "#0000FF", "#006400"))+ggtitle("mtDNA coverage for CI treated #cells")
#geni_mt_id=row.names(data_raw)[row.names(data_raw)%in%as.vector(gene.name.id[grep("mt-Rnr",names(gene#.name.id))])]
#sum_mt=apply(data_raw[geni_mt_id,row.names(classe_epi)],2,sum)
#cluster_condition=factor(cluster_condition,levels=c("Win Epi CI","Intermediate Epi CI","Los Epi CI"))
#data_plot=data.frame(cluster_condition,sum_mt)
#p<-ggplot(data_plot, aes(x=cluster_condition, y=sum_mt, fill=cluster_condition)) +
#geom_boxplot()+labs(fill="Cluster")+ylab("mtDNA content")+xlab("Cluster ")
#p+scale_fill_manual(values=c("#DD6400", "#0000FF", "#006400"))+ggtitle("mtDNA content from only #mt-Rnr1 and mt-Rnr2")
Extended Figure 1 E
data_tpm_log=log10(data_tpm+1)
# mr of repeats for bootsraping
nrRepeats = 100
# percent of gHVGs to use
percent = 60
# nr of HVGs left
sampleSize <- round(length(HVGC)*percent/100)
myBootstrap <- function(indexVector){
idxMatrix <- matrix(rep(NA, sampleSize*nrRepeats)
,nrow = nrRepeats)
for (i in 1:nrow(idxMatrix)){
mySample <- sample(indexVector
,size = sampleSize)
idxMatrix[i,] <- mySample
}
idxMatrix
}
# indices of all HVGs
HVGidx <- which(rownames(data_tpm_log) %in% HVGC)
bootstrapIdxMatrix <- myBootstrap(HVGidx)
# function to calculated clusters
claculateClusters <- function(geneIndices,deepSplit){
chosen.exprs <- (data_tpm_log)[geneIndices,]
# Calculate distance matrix between cells
dissimilarity <- sqrt((1-cor(chosen.exprs, method = "spearman"))/2)
my.dist <- as.dist(dissimilarity)
my.tree <- hclust(my.dist
, method="ward.D2" # method = "ward.D2" was adopted form the reference https://f1000research.com/articles/5-2122/v1 NOTE: next time take method "average"
# some good reference about the different methods is here: https://cran.r-project.org/web/packages/dendextend/vignettes/Cluster_Analysis.html#the-3-clusters-from-the-complete-method-vs-the-real-species-category
)
# Cut the treee to identify clusters
my.clusters <- cutreeHybrid(my.tree
,distM=as.matrix(my.dist)
,deepSplit = deepSplit
)$label
#names(my.clusters) <- clean_annotation$Run_Lane
# calculate the cluster statistics
clusterStatistics <- cluster.stats(d = my.dist
,clustering = my.clusters)
result <- list("PearsonGamma" = clusterStatistics$pearsongamma
,"avgSilhouetteWidth" = clusterStatistics$avg.silwidth)
result
}
# dummy vector of Pearson
# Cluster the data
# Make hierarchical clustering
# dummy vector of Pearson Gammas for 5 deepSplit values and 100 bootstrap repeats
PearsonGammaVector <- matrix(rep(NA, nrRepeats * 5)
,ncol = 5)
colnames(PearsonGammaVector) <- c(0,1,2,3,4)
# dummy vector of average silhouette width for 5 deepSplit values and 100 bootstrap repeats
avgSilhouetteWidth <- matrix(rep(NA, nrRepeats * 5)
,ncol = 5)
colnames(avgSilhouetteWidth) <- c(0,1,2,3,4)
for(j in 0:4){
for (i in 1:nrow(bootstrapIdxMatrix)){
print(paste("i is", i, ", j is" ,j)) # just to keep track of the iterations???
# take indices of HVGs from the bootstrapIdxMatrix of the i-th iteration
indices <- bootstrapIdxMatrix[i,]
# calculate the cluster statistics with these indices and deepSplit == j
myResults <- claculateClusters(indices,j)
# save the Pearsom gammma results
PearsonGammaVector[i,j+1] <- myResults$PearsonGamma
# save the average silhouette width results
avgSilhouetteWidth[i,j+1] <- myResults$avgSilhouetteWidth
}
}
## [1] "i is 1 , j is 0"
## ..cutHeight not given, setting it to 4.22 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 2 , j is 0"
## ..cutHeight not given, setting it to 4.2 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 3 , j is 0"
## ..cutHeight not given, setting it to 4.38 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 4 , j is 0"
## ..cutHeight not given, setting it to 4.19 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 5 , j is 0"
## ..cutHeight not given, setting it to 4.11 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 6 , j is 0"
## ..cutHeight not given, setting it to 4.34 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 7 , j is 0"
## ..cutHeight not given, setting it to 4.4 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 8 , j is 0"
## ..cutHeight not given, setting it to 4.12 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 9 , j is 0"
## ..cutHeight not given, setting it to 4.24 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 10 , j is 0"
## ..cutHeight not given, setting it to 4.27 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 11 , j is 0"
## ..cutHeight not given, setting it to 4.34 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 12 , j is 0"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 13 , j is 0"
## ..cutHeight not given, setting it to 4.2 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 14 , j is 0"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 15 , j is 0"
## ..cutHeight not given, setting it to 4.33 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 16 , j is 0"
## ..cutHeight not given, setting it to 4.08 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 17 , j is 0"
## ..cutHeight not given, setting it to 4.37 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 18 , j is 0"
## ..cutHeight not given, setting it to 4.21 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 19 , j is 0"
## ..cutHeight not given, setting it to 4.21 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 20 , j is 0"
## ..cutHeight not given, setting it to 4.27 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 21 , j is 0"
## ..cutHeight not given, setting it to 4.05 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 22 , j is 0"
## ..cutHeight not given, setting it to 4.37 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 23 , j is 0"
## ..cutHeight not given, setting it to 4.16 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 24 , j is 0"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 25 , j is 0"
## ..cutHeight not given, setting it to 4.22 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 26 , j is 0"
## ..cutHeight not given, setting it to 4.25 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 27 , j is 0"
## ..cutHeight not given, setting it to 4.29 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 28 , j is 0"
## ..cutHeight not given, setting it to 4.17 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 29 , j is 0"
## ..cutHeight not given, setting it to 4.21 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 30 , j is 0"
## ..cutHeight not given, setting it to 4.38 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 31 , j is 0"
## ..cutHeight not given, setting it to 4.38 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 32 , j is 0"
## ..cutHeight not given, setting it to 4.22 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 33 , j is 0"
## ..cutHeight not given, setting it to 4.39 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 34 , j is 0"
## ..cutHeight not given, setting it to 4.33 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 35 , j is 0"
## ..cutHeight not given, setting it to 4.45 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 36 , j is 0"
## ..cutHeight not given, setting it to 4.23 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 37 , j is 0"
## ..cutHeight not given, setting it to 4.25 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 38 , j is 0"
## ..cutHeight not given, setting it to 4.23 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 39 , j is 0"
## ..cutHeight not given, setting it to 4.34 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 40 , j is 0"
## ..cutHeight not given, setting it to 4.12 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 41 , j is 0"
## ..cutHeight not given, setting it to 4.15 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 42 , j is 0"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 43 , j is 0"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 44 , j is 0"
## ..cutHeight not given, setting it to 4.27 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 45 , j is 0"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 46 , j is 0"
## ..cutHeight not given, setting it to 4.44 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 47 , j is 0"
## ..cutHeight not given, setting it to 4.25 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 48 , j is 0"
## ..cutHeight not given, setting it to 4.21 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 49 , j is 0"
## ..cutHeight not given, setting it to 4.23 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 50 , j is 0"
## ..cutHeight not given, setting it to 4.22 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 51 , j is 0"
## ..cutHeight not given, setting it to 4.34 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 52 , j is 0"
## ..cutHeight not given, setting it to 4.16 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 53 , j is 0"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 54 , j is 0"
## ..cutHeight not given, setting it to 4.16 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 55 , j is 0"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 56 , j is 0"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 57 , j is 0"
## ..cutHeight not given, setting it to 4.26 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 58 , j is 0"
## ..cutHeight not given, setting it to 4.09 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 59 , j is 0"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 60 , j is 0"
## ..cutHeight not given, setting it to 4.43 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 61 , j is 0"
## ..cutHeight not given, setting it to 4.39 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 62 , j is 0"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 63 , j is 0"
## ..cutHeight not given, setting it to 4.28 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 64 , j is 0"
## ..cutHeight not given, setting it to 4.17 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 65 , j is 0"
## ..cutHeight not given, setting it to 4.18 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 66 , j is 0"
## ..cutHeight not given, setting it to 4.29 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 67 , j is 0"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 68 , j is 0"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 69 , j is 0"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 70 , j is 0"
## ..cutHeight not given, setting it to 4.26 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 71 , j is 0"
## ..cutHeight not given, setting it to 4.21 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 72 , j is 0"
## ..cutHeight not given, setting it to 4.27 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 73 , j is 0"
## ..cutHeight not given, setting it to 4.25 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 74 , j is 0"
## ..cutHeight not given, setting it to 4.1 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 75 , j is 0"
## ..cutHeight not given, setting it to 4.29 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 76 , j is 0"
## ..cutHeight not given, setting it to 4.28 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 77 , j is 0"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 78 , j is 0"
## ..cutHeight not given, setting it to 4.52 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 79 , j is 0"
## ..cutHeight not given, setting it to 4.33 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 80 , j is 0"
## ..cutHeight not given, setting it to 4.46 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 81 , j is 0"
## ..cutHeight not given, setting it to 4.42 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 82 , j is 0"
## ..cutHeight not given, setting it to 4.25 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 83 , j is 0"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 84 , j is 0"
## ..cutHeight not given, setting it to 4.33 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 85 , j is 0"
## ..cutHeight not given, setting it to 4.19 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 86 , j is 0"
## ..cutHeight not given, setting it to 4.39 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 87 , j is 0"
## ..cutHeight not given, setting it to 4.4 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 88 , j is 0"
## ..cutHeight not given, setting it to 4.39 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 89 , j is 0"
## ..cutHeight not given, setting it to 4.34 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 90 , j is 0"
## ..cutHeight not given, setting it to 4.26 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 91 , j is 0"
## ..cutHeight not given, setting it to 4.5 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 92 , j is 0"
## ..cutHeight not given, setting it to 4.4 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 93 , j is 0"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 94 , j is 0"
## ..cutHeight not given, setting it to 4.39 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 95 , j is 0"
## ..cutHeight not given, setting it to 4.29 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 96 , j is 0"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 97 , j is 0"
## ..cutHeight not given, setting it to 4.29 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 98 , j is 0"
## ..cutHeight not given, setting it to 4.4 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 99 , j is 0"
## ..cutHeight not given, setting it to 4.38 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 100 , j is 0"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 1 , j is 1"
## ..cutHeight not given, setting it to 4.22 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 2 , j is 1"
## ..cutHeight not given, setting it to 4.2 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 3 , j is 1"
## ..cutHeight not given, setting it to 4.38 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 4 , j is 1"
## ..cutHeight not given, setting it to 4.19 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 5 , j is 1"
## ..cutHeight not given, setting it to 4.11 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 6 , j is 1"
## ..cutHeight not given, setting it to 4.34 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 7 , j is 1"
## ..cutHeight not given, setting it to 4.4 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 8 , j is 1"
## ..cutHeight not given, setting it to 4.12 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 9 , j is 1"
## ..cutHeight not given, setting it to 4.24 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 10 , j is 1"
## ..cutHeight not given, setting it to 4.27 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 11 , j is 1"
## ..cutHeight not given, setting it to 4.34 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 12 , j is 1"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 13 , j is 1"
## ..cutHeight not given, setting it to 4.2 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 14 , j is 1"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 15 , j is 1"
## ..cutHeight not given, setting it to 4.33 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 16 , j is 1"
## ..cutHeight not given, setting it to 4.08 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 17 , j is 1"
## ..cutHeight not given, setting it to 4.37 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 18 , j is 1"
## ..cutHeight not given, setting it to 4.21 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 19 , j is 1"
## ..cutHeight not given, setting it to 4.21 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 20 , j is 1"
## ..cutHeight not given, setting it to 4.27 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 21 , j is 1"
## ..cutHeight not given, setting it to 4.05 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 22 , j is 1"
## ..cutHeight not given, setting it to 4.37 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 23 , j is 1"
## ..cutHeight not given, setting it to 4.16 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 24 , j is 1"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 25 , j is 1"
## ..cutHeight not given, setting it to 4.22 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 26 , j is 1"
## ..cutHeight not given, setting it to 4.25 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 27 , j is 1"
## ..cutHeight not given, setting it to 4.29 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 28 , j is 1"
## ..cutHeight not given, setting it to 4.17 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 29 , j is 1"
## ..cutHeight not given, setting it to 4.21 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 30 , j is 1"
## ..cutHeight not given, setting it to 4.38 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 31 , j is 1"
## ..cutHeight not given, setting it to 4.38 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 32 , j is 1"
## ..cutHeight not given, setting it to 4.22 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 33 , j is 1"
## ..cutHeight not given, setting it to 4.39 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 34 , j is 1"
## ..cutHeight not given, setting it to 4.33 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 35 , j is 1"
## ..cutHeight not given, setting it to 4.45 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 36 , j is 1"
## ..cutHeight not given, setting it to 4.23 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 37 , j is 1"
## ..cutHeight not given, setting it to 4.25 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 38 , j is 1"
## ..cutHeight not given, setting it to 4.23 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 39 , j is 1"
## ..cutHeight not given, setting it to 4.34 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 40 , j is 1"
## ..cutHeight not given, setting it to 4.12 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 41 , j is 1"
## ..cutHeight not given, setting it to 4.15 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 42 , j is 1"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 43 , j is 1"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 44 , j is 1"
## ..cutHeight not given, setting it to 4.27 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 45 , j is 1"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 46 , j is 1"
## ..cutHeight not given, setting it to 4.44 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 47 , j is 1"
## ..cutHeight not given, setting it to 4.25 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 48 , j is 1"
## ..cutHeight not given, setting it to 4.21 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 49 , j is 1"
## ..cutHeight not given, setting it to 4.23 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 50 , j is 1"
## ..cutHeight not given, setting it to 4.22 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 51 , j is 1"
## ..cutHeight not given, setting it to 4.34 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 52 , j is 1"
## ..cutHeight not given, setting it to 4.16 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 53 , j is 1"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 54 , j is 1"
## ..cutHeight not given, setting it to 4.16 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 55 , j is 1"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 56 , j is 1"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 57 , j is 1"
## ..cutHeight not given, setting it to 4.26 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 58 , j is 1"
## ..cutHeight not given, setting it to 4.09 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 59 , j is 1"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 60 , j is 1"
## ..cutHeight not given, setting it to 4.43 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 61 , j is 1"
## ..cutHeight not given, setting it to 4.39 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 62 , j is 1"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 63 , j is 1"
## ..cutHeight not given, setting it to 4.28 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 64 , j is 1"
## ..cutHeight not given, setting it to 4.17 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 65 , j is 1"
## ..cutHeight not given, setting it to 4.18 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 66 , j is 1"
## ..cutHeight not given, setting it to 4.29 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 67 , j is 1"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 68 , j is 1"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 69 , j is 1"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 70 , j is 1"
## ..cutHeight not given, setting it to 4.26 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 71 , j is 1"
## ..cutHeight not given, setting it to 4.21 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 72 , j is 1"
## ..cutHeight not given, setting it to 4.27 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 73 , j is 1"
## ..cutHeight not given, setting it to 4.25 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 74 , j is 1"
## ..cutHeight not given, setting it to 4.1 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 75 , j is 1"
## ..cutHeight not given, setting it to 4.29 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 76 , j is 1"
## ..cutHeight not given, setting it to 4.28 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 77 , j is 1"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 78 , j is 1"
## ..cutHeight not given, setting it to 4.52 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 79 , j is 1"
## ..cutHeight not given, setting it to 4.33 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 80 , j is 1"
## ..cutHeight not given, setting it to 4.46 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 81 , j is 1"
## ..cutHeight not given, setting it to 4.42 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 82 , j is 1"
## ..cutHeight not given, setting it to 4.25 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 83 , j is 1"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 84 , j is 1"
## ..cutHeight not given, setting it to 4.33 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 85 , j is 1"
## ..cutHeight not given, setting it to 4.19 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 86 , j is 1"
## ..cutHeight not given, setting it to 4.39 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 87 , j is 1"
## ..cutHeight not given, setting it to 4.4 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 88 , j is 1"
## ..cutHeight not given, setting it to 4.39 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 89 , j is 1"
## ..cutHeight not given, setting it to 4.34 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 90 , j is 1"
## ..cutHeight not given, setting it to 4.26 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 91 , j is 1"
## ..cutHeight not given, setting it to 4.5 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 92 , j is 1"
## ..cutHeight not given, setting it to 4.4 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 93 , j is 1"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 94 , j is 1"
## ..cutHeight not given, setting it to 4.39 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 95 , j is 1"
## ..cutHeight not given, setting it to 4.29 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 96 , j is 1"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 97 , j is 1"
## ..cutHeight not given, setting it to 4.29 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 98 , j is 1"
## ..cutHeight not given, setting it to 4.4 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 99 , j is 1"
## ..cutHeight not given, setting it to 4.38 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 100 , j is 1"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 1 , j is 2"
## ..cutHeight not given, setting it to 4.22 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 2 , j is 2"
## ..cutHeight not given, setting it to 4.2 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 3 , j is 2"
## ..cutHeight not given, setting it to 4.38 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 4 , j is 2"
## ..cutHeight not given, setting it to 4.19 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 5 , j is 2"
## ..cutHeight not given, setting it to 4.11 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 6 , j is 2"
## ..cutHeight not given, setting it to 4.34 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 7 , j is 2"
## ..cutHeight not given, setting it to 4.4 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 8 , j is 2"
## ..cutHeight not given, setting it to 4.12 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 9 , j is 2"
## ..cutHeight not given, setting it to 4.24 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 10 , j is 2"
## ..cutHeight not given, setting it to 4.27 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 11 , j is 2"
## ..cutHeight not given, setting it to 4.34 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 12 , j is 2"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 13 , j is 2"
## ..cutHeight not given, setting it to 4.2 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 14 , j is 2"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 15 , j is 2"
## ..cutHeight not given, setting it to 4.33 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 16 , j is 2"
## ..cutHeight not given, setting it to 4.08 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 17 , j is 2"
## ..cutHeight not given, setting it to 4.37 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 18 , j is 2"
## ..cutHeight not given, setting it to 4.21 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 19 , j is 2"
## ..cutHeight not given, setting it to 4.21 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 20 , j is 2"
## ..cutHeight not given, setting it to 4.27 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 21 , j is 2"
## ..cutHeight not given, setting it to 4.05 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 22 , j is 2"
## ..cutHeight not given, setting it to 4.37 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 23 , j is 2"
## ..cutHeight not given, setting it to 4.16 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 24 , j is 2"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 25 , j is 2"
## ..cutHeight not given, setting it to 4.22 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 26 , j is 2"
## ..cutHeight not given, setting it to 4.25 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 27 , j is 2"
## ..cutHeight not given, setting it to 4.29 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 28 , j is 2"
## ..cutHeight not given, setting it to 4.17 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 29 , j is 2"
## ..cutHeight not given, setting it to 4.21 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 30 , j is 2"
## ..cutHeight not given, setting it to 4.38 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 31 , j is 2"
## ..cutHeight not given, setting it to 4.38 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 32 , j is 2"
## ..cutHeight not given, setting it to 4.22 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 33 , j is 2"
## ..cutHeight not given, setting it to 4.39 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 34 , j is 2"
## ..cutHeight not given, setting it to 4.33 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 35 , j is 2"
## ..cutHeight not given, setting it to 4.45 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 36 , j is 2"
## ..cutHeight not given, setting it to 4.23 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 37 , j is 2"
## ..cutHeight not given, setting it to 4.25 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 38 , j is 2"
## ..cutHeight not given, setting it to 4.23 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 39 , j is 2"
## ..cutHeight not given, setting it to 4.34 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 40 , j is 2"
## ..cutHeight not given, setting it to 4.12 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 41 , j is 2"
## ..cutHeight not given, setting it to 4.15 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 42 , j is 2"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 43 , j is 2"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 44 , j is 2"
## ..cutHeight not given, setting it to 4.27 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 45 , j is 2"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 46 , j is 2"
## ..cutHeight not given, setting it to 4.44 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 47 , j is 2"
## ..cutHeight not given, setting it to 4.25 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 48 , j is 2"
## ..cutHeight not given, setting it to 4.21 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 49 , j is 2"
## ..cutHeight not given, setting it to 4.23 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 50 , j is 2"
## ..cutHeight not given, setting it to 4.22 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 51 , j is 2"
## ..cutHeight not given, setting it to 4.34 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 52 , j is 2"
## ..cutHeight not given, setting it to 4.16 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 53 , j is 2"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 54 , j is 2"
## ..cutHeight not given, setting it to 4.16 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 55 , j is 2"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 56 , j is 2"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 57 , j is 2"
## ..cutHeight not given, setting it to 4.26 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 58 , j is 2"
## ..cutHeight not given, setting it to 4.09 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 59 , j is 2"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 60 , j is 2"
## ..cutHeight not given, setting it to 4.43 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 61 , j is 2"
## ..cutHeight not given, setting it to 4.39 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 62 , j is 2"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 63 , j is 2"
## ..cutHeight not given, setting it to 4.28 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 64 , j is 2"
## ..cutHeight not given, setting it to 4.17 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 65 , j is 2"
## ..cutHeight not given, setting it to 4.18 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 66 , j is 2"
## ..cutHeight not given, setting it to 4.29 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 67 , j is 2"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 68 , j is 2"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 69 , j is 2"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 70 , j is 2"
## ..cutHeight not given, setting it to 4.26 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 71 , j is 2"
## ..cutHeight not given, setting it to 4.21 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 72 , j is 2"
## ..cutHeight not given, setting it to 4.27 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 73 , j is 2"
## ..cutHeight not given, setting it to 4.25 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 74 , j is 2"
## ..cutHeight not given, setting it to 4.1 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 75 , j is 2"
## ..cutHeight not given, setting it to 4.29 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 76 , j is 2"
## ..cutHeight not given, setting it to 4.28 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 77 , j is 2"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 78 , j is 2"
## ..cutHeight not given, setting it to 4.52 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 79 , j is 2"
## ..cutHeight not given, setting it to 4.33 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 80 , j is 2"
## ..cutHeight not given, setting it to 4.46 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 81 , j is 2"
## ..cutHeight not given, setting it to 4.42 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 82 , j is 2"
## ..cutHeight not given, setting it to 4.25 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 83 , j is 2"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 84 , j is 2"
## ..cutHeight not given, setting it to 4.33 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 85 , j is 2"
## ..cutHeight not given, setting it to 4.19 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 86 , j is 2"
## ..cutHeight not given, setting it to 4.39 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 87 , j is 2"
## ..cutHeight not given, setting it to 4.4 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 88 , j is 2"
## ..cutHeight not given, setting it to 4.39 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 89 , j is 2"
## ..cutHeight not given, setting it to 4.34 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 90 , j is 2"
## ..cutHeight not given, setting it to 4.26 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 91 , j is 2"
## ..cutHeight not given, setting it to 4.5 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 92 , j is 2"
## ..cutHeight not given, setting it to 4.4 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 93 , j is 2"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 94 , j is 2"
## ..cutHeight not given, setting it to 4.39 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 95 , j is 2"
## ..cutHeight not given, setting it to 4.29 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 96 , j is 2"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 97 , j is 2"
## ..cutHeight not given, setting it to 4.29 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 98 , j is 2"
## ..cutHeight not given, setting it to 4.4 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 99 , j is 2"
## ..cutHeight not given, setting it to 4.38 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 100 , j is 2"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 1 , j is 3"
## ..cutHeight not given, setting it to 4.22 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 2 , j is 3"
## ..cutHeight not given, setting it to 4.2 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 3 , j is 3"
## ..cutHeight not given, setting it to 4.38 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 4 , j is 3"
## ..cutHeight not given, setting it to 4.19 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 5 , j is 3"
## ..cutHeight not given, setting it to 4.11 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 6 , j is 3"
## ..cutHeight not given, setting it to 4.34 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 7 , j is 3"
## ..cutHeight not given, setting it to 4.4 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 8 , j is 3"
## ..cutHeight not given, setting it to 4.12 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 9 , j is 3"
## ..cutHeight not given, setting it to 4.24 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 10 , j is 3"
## ..cutHeight not given, setting it to 4.27 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 11 , j is 3"
## ..cutHeight not given, setting it to 4.34 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 12 , j is 3"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 13 , j is 3"
## ..cutHeight not given, setting it to 4.2 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 14 , j is 3"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 15 , j is 3"
## ..cutHeight not given, setting it to 4.33 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 16 , j is 3"
## ..cutHeight not given, setting it to 4.08 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 17 , j is 3"
## ..cutHeight not given, setting it to 4.37 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 18 , j is 3"
## ..cutHeight not given, setting it to 4.21 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 19 , j is 3"
## ..cutHeight not given, setting it to 4.21 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 20 , j is 3"
## ..cutHeight not given, setting it to 4.27 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 21 , j is 3"
## ..cutHeight not given, setting it to 4.05 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 22 , j is 3"
## ..cutHeight not given, setting it to 4.37 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 23 , j is 3"
## ..cutHeight not given, setting it to 4.16 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 24 , j is 3"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 25 , j is 3"
## ..cutHeight not given, setting it to 4.22 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 26 , j is 3"
## ..cutHeight not given, setting it to 4.25 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 27 , j is 3"
## ..cutHeight not given, setting it to 4.29 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 28 , j is 3"
## ..cutHeight not given, setting it to 4.17 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 29 , j is 3"
## ..cutHeight not given, setting it to 4.21 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 30 , j is 3"
## ..cutHeight not given, setting it to 4.38 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 31 , j is 3"
## ..cutHeight not given, setting it to 4.38 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 32 , j is 3"
## ..cutHeight not given, setting it to 4.22 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 33 , j is 3"
## ..cutHeight not given, setting it to 4.39 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 34 , j is 3"
## ..cutHeight not given, setting it to 4.33 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 35 , j is 3"
## ..cutHeight not given, setting it to 4.45 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 36 , j is 3"
## ..cutHeight not given, setting it to 4.23 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 37 , j is 3"
## ..cutHeight not given, setting it to 4.25 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 38 , j is 3"
## ..cutHeight not given, setting it to 4.23 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 39 , j is 3"
## ..cutHeight not given, setting it to 4.34 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 40 , j is 3"
## ..cutHeight not given, setting it to 4.12 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 41 , j is 3"
## ..cutHeight not given, setting it to 4.15 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 42 , j is 3"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 43 , j is 3"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 44 , j is 3"
## ..cutHeight not given, setting it to 4.27 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 45 , j is 3"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 46 , j is 3"
## ..cutHeight not given, setting it to 4.44 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 47 , j is 3"
## ..cutHeight not given, setting it to 4.25 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 48 , j is 3"
## ..cutHeight not given, setting it to 4.21 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 49 , j is 3"
## ..cutHeight not given, setting it to 4.23 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 50 , j is 3"
## ..cutHeight not given, setting it to 4.22 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 51 , j is 3"
## ..cutHeight not given, setting it to 4.34 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 52 , j is 3"
## ..cutHeight not given, setting it to 4.16 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 53 , j is 3"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 54 , j is 3"
## ..cutHeight not given, setting it to 4.16 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 55 , j is 3"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 56 , j is 3"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 57 , j is 3"
## ..cutHeight not given, setting it to 4.26 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 58 , j is 3"
## ..cutHeight not given, setting it to 4.09 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 59 , j is 3"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 60 , j is 3"
## ..cutHeight not given, setting it to 4.43 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 61 , j is 3"
## ..cutHeight not given, setting it to 4.39 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 62 , j is 3"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 63 , j is 3"
## ..cutHeight not given, setting it to 4.28 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 64 , j is 3"
## ..cutHeight not given, setting it to 4.17 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 65 , j is 3"
## ..cutHeight not given, setting it to 4.18 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 66 , j is 3"
## ..cutHeight not given, setting it to 4.29 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 67 , j is 3"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 68 , j is 3"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 69 , j is 3"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 70 , j is 3"
## ..cutHeight not given, setting it to 4.26 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 71 , j is 3"
## ..cutHeight not given, setting it to 4.21 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 72 , j is 3"
## ..cutHeight not given, setting it to 4.27 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 73 , j is 3"
## ..cutHeight not given, setting it to 4.25 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 74 , j is 3"
## ..cutHeight not given, setting it to 4.1 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 75 , j is 3"
## ..cutHeight not given, setting it to 4.29 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 76 , j is 3"
## ..cutHeight not given, setting it to 4.28 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 77 , j is 3"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 78 , j is 3"
## ..cutHeight not given, setting it to 4.52 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 79 , j is 3"
## ..cutHeight not given, setting it to 4.33 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 80 , j is 3"
## ..cutHeight not given, setting it to 4.46 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 81 , j is 3"
## ..cutHeight not given, setting it to 4.42 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 82 , j is 3"
## ..cutHeight not given, setting it to 4.25 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 83 , j is 3"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 84 , j is 3"
## ..cutHeight not given, setting it to 4.33 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 85 , j is 3"
## ..cutHeight not given, setting it to 4.19 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 86 , j is 3"
## ..cutHeight not given, setting it to 4.39 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 87 , j is 3"
## ..cutHeight not given, setting it to 4.4 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 88 , j is 3"
## ..cutHeight not given, setting it to 4.39 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 89 , j is 3"
## ..cutHeight not given, setting it to 4.34 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 90 , j is 3"
## ..cutHeight not given, setting it to 4.26 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 91 , j is 3"
## ..cutHeight not given, setting it to 4.5 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 92 , j is 3"
## ..cutHeight not given, setting it to 4.4 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 93 , j is 3"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 94 , j is 3"
## ..cutHeight not given, setting it to 4.39 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 95 , j is 3"
## ..cutHeight not given, setting it to 4.29 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 96 , j is 3"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 97 , j is 3"
## ..cutHeight not given, setting it to 4.29 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 98 , j is 3"
## ..cutHeight not given, setting it to 4.4 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 99 , j is 3"
## ..cutHeight not given, setting it to 4.38 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 100 , j is 3"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 1 , j is 4"
## ..cutHeight not given, setting it to 4.22 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 2 , j is 4"
## ..cutHeight not given, setting it to 4.2 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 3 , j is 4"
## ..cutHeight not given, setting it to 4.38 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 4 , j is 4"
## ..cutHeight not given, setting it to 4.19 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 5 , j is 4"
## ..cutHeight not given, setting it to 4.11 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 6 , j is 4"
## ..cutHeight not given, setting it to 4.34 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 7 , j is 4"
## ..cutHeight not given, setting it to 4.4 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 8 , j is 4"
## ..cutHeight not given, setting it to 4.12 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 9 , j is 4"
## ..cutHeight not given, setting it to 4.24 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 10 , j is 4"
## ..cutHeight not given, setting it to 4.27 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 11 , j is 4"
## ..cutHeight not given, setting it to 4.34 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 12 , j is 4"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 13 , j is 4"
## ..cutHeight not given, setting it to 4.2 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 14 , j is 4"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 15 , j is 4"
## ..cutHeight not given, setting it to 4.33 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 16 , j is 4"
## ..cutHeight not given, setting it to 4.08 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 17 , j is 4"
## ..cutHeight not given, setting it to 4.37 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 18 , j is 4"
## ..cutHeight not given, setting it to 4.21 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 19 , j is 4"
## ..cutHeight not given, setting it to 4.21 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 20 , j is 4"
## ..cutHeight not given, setting it to 4.27 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 21 , j is 4"
## ..cutHeight not given, setting it to 4.05 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 22 , j is 4"
## ..cutHeight not given, setting it to 4.37 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 23 , j is 4"
## ..cutHeight not given, setting it to 4.16 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 24 , j is 4"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 25 , j is 4"
## ..cutHeight not given, setting it to 4.22 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 26 , j is 4"
## ..cutHeight not given, setting it to 4.25 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 27 , j is 4"
## ..cutHeight not given, setting it to 4.29 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 28 , j is 4"
## ..cutHeight not given, setting it to 4.17 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 29 , j is 4"
## ..cutHeight not given, setting it to 4.21 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 30 , j is 4"
## ..cutHeight not given, setting it to 4.38 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 31 , j is 4"
## ..cutHeight not given, setting it to 4.38 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 32 , j is 4"
## ..cutHeight not given, setting it to 4.22 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 33 , j is 4"
## ..cutHeight not given, setting it to 4.39 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 34 , j is 4"
## ..cutHeight not given, setting it to 4.33 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 35 , j is 4"
## ..cutHeight not given, setting it to 4.45 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 36 , j is 4"
## ..cutHeight not given, setting it to 4.23 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 37 , j is 4"
## ..cutHeight not given, setting it to 4.25 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 38 , j is 4"
## ..cutHeight not given, setting it to 4.23 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 39 , j is 4"
## ..cutHeight not given, setting it to 4.34 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 40 , j is 4"
## ..cutHeight not given, setting it to 4.12 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 41 , j is 4"
## ..cutHeight not given, setting it to 4.15 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 42 , j is 4"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 43 , j is 4"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 44 , j is 4"
## ..cutHeight not given, setting it to 4.27 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 45 , j is 4"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 46 , j is 4"
## ..cutHeight not given, setting it to 4.44 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 47 , j is 4"
## ..cutHeight not given, setting it to 4.25 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 48 , j is 4"
## ..cutHeight not given, setting it to 4.21 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 49 , j is 4"
## ..cutHeight not given, setting it to 4.23 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 50 , j is 4"
## ..cutHeight not given, setting it to 4.22 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 51 , j is 4"
## ..cutHeight not given, setting it to 4.34 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 52 , j is 4"
## ..cutHeight not given, setting it to 4.16 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 53 , j is 4"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 54 , j is 4"
## ..cutHeight not given, setting it to 4.16 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 55 , j is 4"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 56 , j is 4"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 57 , j is 4"
## ..cutHeight not given, setting it to 4.26 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 58 , j is 4"
## ..cutHeight not given, setting it to 4.09 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 59 , j is 4"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 60 , j is 4"
## ..cutHeight not given, setting it to 4.43 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 61 , j is 4"
## ..cutHeight not given, setting it to 4.39 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 62 , j is 4"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 63 , j is 4"
## ..cutHeight not given, setting it to 4.28 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 64 , j is 4"
## ..cutHeight not given, setting it to 4.17 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 65 , j is 4"
## ..cutHeight not given, setting it to 4.18 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 66 , j is 4"
## ..cutHeight not given, setting it to 4.29 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 67 , j is 4"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 68 , j is 4"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 69 , j is 4"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 70 , j is 4"
## ..cutHeight not given, setting it to 4.26 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 71 , j is 4"
## ..cutHeight not given, setting it to 4.21 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 72 , j is 4"
## ..cutHeight not given, setting it to 4.27 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 73 , j is 4"
## ..cutHeight not given, setting it to 4.25 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 74 , j is 4"
## ..cutHeight not given, setting it to 4.1 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 75 , j is 4"
## ..cutHeight not given, setting it to 4.29 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 76 , j is 4"
## ..cutHeight not given, setting it to 4.28 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 77 , j is 4"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 78 , j is 4"
## ..cutHeight not given, setting it to 4.52 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 79 , j is 4"
## ..cutHeight not given, setting it to 4.33 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 80 , j is 4"
## ..cutHeight not given, setting it to 4.46 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 81 , j is 4"
## ..cutHeight not given, setting it to 4.42 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 82 , j is 4"
## ..cutHeight not given, setting it to 4.25 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 83 , j is 4"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 84 , j is 4"
## ..cutHeight not given, setting it to 4.33 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 85 , j is 4"
## ..cutHeight not given, setting it to 4.19 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 86 , j is 4"
## ..cutHeight not given, setting it to 4.39 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 87 , j is 4"
## ..cutHeight not given, setting it to 4.4 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 88 , j is 4"
## ..cutHeight not given, setting it to 4.39 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 89 , j is 4"
## ..cutHeight not given, setting it to 4.34 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 90 , j is 4"
## ..cutHeight not given, setting it to 4.26 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 91 , j is 4"
## ..cutHeight not given, setting it to 4.5 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 92 , j is 4"
## ..cutHeight not given, setting it to 4.4 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 93 , j is 4"
## ..cutHeight not given, setting it to 4.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 94 , j is 4"
## ..cutHeight not given, setting it to 4.39 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 95 , j is 4"
## ..cutHeight not given, setting it to 4.29 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 96 , j is 4"
## ..cutHeight not given, setting it to 4.32 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 97 , j is 4"
## ..cutHeight not given, setting it to 4.29 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 98 , j is 4"
## ..cutHeight not given, setting it to 4.4 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 99 , j is 4"
## ..cutHeight not given, setting it to 4.38 ===> 99% of the (truncated) height range in dendro.
## ..done.
## [1] "i is 100 , j is 4"
## ..cutHeight not given, setting it to 4.36 ===> 99% of the (truncated) height range in dendro.
## ..done.
pearsonGammaBoxplot <- boxplot(PearsonGammaVector
,use.cols = TRUE
,main = "Pearson Gamma"
,xlab = "deepSplit"
,ylab = "Pearson Gamma"
)

pearsonGammaBoxplot
## $stats
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.6453768 0.6746437 0.4780571 0.2711807 0.2292184
## [2,] 0.6725865 0.7028253 0.5210694 0.2972013 0.2571013
## [3,] 0.6848418 0.7118417 0.5474299 0.3108514 0.2678744
## [4,] 0.6960932 0.7218877 0.5822492 0.3233165 0.2768857
## [5,] 0.7213358 0.7387877 0.6646111 0.3606915 0.2997035
##
## $n
## [1] 100 100 100 100 100
##
## $conf
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.6811277 0.7088298 0.5377635 0.3067252 0.2647485
## [2,] 0.6885558 0.7148535 0.5570963 0.3149776 0.2710004
##
## $out
## [1] 0.6725404 0.7065068 0.7074678 0.7095835 0.7153105 0.3627821 0.2577033
## [8] 0.3659581 0.3913694 0.3103310 0.2264069
##
## $group
## [1] 2 3 3 3 3 4 4 4 4 5 5
##
## $names
## [1] "0" "1" "2" "3" "4"
avgSilWidth <- boxplot(avgSilhouetteWidth
,use.cols = TRUE
,main = "Average Silhouette Width"
,xlab = "deepSplit"
,ylab = "Average Silhouette Width"
)

avgSilWidth
## $stats
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.1150810 0.1167279 0.05311246 0.007183185 0.002651091
## [2,] 0.1207772 0.1218660 0.05973173 0.011508356 0.008094263
## [3,] 0.1245282 0.1252157 0.06217317 0.013632014 0.010221661
## [4,] 0.1270027 0.1271003 0.06463196 0.015648046 0.011810760
## [5,] 0.1359728 0.1333532 0.07030367 0.020342627 0.016083051
##
## $n
## [1] 100 100 100 100 100
##
## $conf
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.1235445 0.1243887 0.06139894 0.01297794 0.009634454
## [2,] 0.1255118 0.1260427 0.06294741 0.01428609 0.010808867
##
## $out
## [1] 0.1104956767 0.1098077069 0.1119631908 0.1178714581 0.1274893688
## [6] 0.0437664888 0.0728402279 0.0437710142 0.0411813026 0.0443882455
## [11] 0.0728568794 0.1032799902 0.0464366720 0.1151405182 0.0419362484
## [16] 0.0802955256 0.0019551496 0.0008981727 0.0024309241
##
## $group
## [1] 1 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 5 5 5
##
## $names
## [1] "0" "1" "2" "3" "4"
Figure 1 B
my.clusters=as.vector(my.clusters)
index_1=grep("1",my.clusters)
index_2=grep("2",my.clusters)
index_3=grep("3",my.clusters)
index_4=grep("4",my.clusters)
index_5=grep("5",my.clusters)
condition_1=Condition[index_1]
condition_1_ci=grep("Cell competition OFF",condition_1)
condition_1_dmso=grep("Cell competition ON",condition_1)
index_1_fine=c(condition_1_ci,condition_1_dmso)
medium_1=condition_1[index_1_fine]
condition_2=Condition[index_2]
condition_2_ci=grep("Cell competition OFF",condition_2)
condition_2_dmso=grep("Cell competition ON",condition_2)
index_2_fine=c(condition_2_ci,condition_2_dmso)
medium_2=condition_2[index_2_fine]
condition_3=Condition[index_3]
condition_3_ci=grep("Cell competition OFF",condition_3)
condition_3_dmso=grep("Cell competition ON",condition_3)
index_3_fine=c(condition_3_ci,condition_3_dmso)
medium_3=condition_3[index_3_fine]
condition_4=Condition[index_4]
condition_4_ci=grep("Cell competition OFF",condition_4)
condition_4_dmso=grep("Cell competition ON",condition_4)
index_4_fine=c(condition_4_ci,condition_4_dmso)
medium_4=condition_4[index_4_fine]
condition_5=Condition[index_5]
condition_5_ci=grep("Cell competition OFF",condition_5)
condition_5_dmso=grep("Cell competition ON",condition_5)
index_5_fine=c(condition_5_ci,condition_5_dmso)
medium_5=condition_5[index_5_fine]
good_1=good[index_1]
good_finale_1=good_1[index_1_fine]
good_2=good[index_2]
good_finale_2=good_2[index_2_fine]
good_3=good[index_3]
good_finale_3=good_3[index_3_fine]
good_4=good[index_4]
good_finale_4=good_4[index_4_fine]
good_5=good[index_5]
good_finale_5=good_5[index_5_fine]
cluster_1=my.clusters[index_1]
cluster_finale_1=cluster_1[index_1_fine]
cluster_2=my.clusters[index_2]
cluster_finale_2=cluster_2[index_2_fine]
cluster_3=my.clusters[index_3]
cluster_finale_3=cluster_3[index_3_fine]
cluster_4=my.clusters[index_4]
cluster_finale_4=cluster_4[index_4_fine]
cluster_5=my.clusters[index_5]
cluster_finale_5=cluster_5[index_5_fine]
cluster_finale=c(cluster_finale_1,cluster_finale_3,cluster_finale_4,cluster_finale_2,cluster_finale_5)
good_finale=c(good_finale_1,good_finale_3,good_finale_4,good_finale_2,good_finale_5)
medium_finale=c(medium_1,medium_3,medium_4,medium_2,medium_5)
geni_icb_exe=c(nome_gene$Gene.stable.ID[which(nome_gene$Gene.name=="Cdx2")],nome_gene$Gene.stable.ID[which(nome_gene$Gene.name=="Bmp4")],nome_gene$Gene.stable.ID[which(nome_gene$Gene.name=="Elf5")])
geni_icb_epi=c(nome_gene$Gene.stable.ID[which(nome_gene$Gene.name=="Gng3")],nome_gene$Gene.stable.ID[which(nome_gene$Gene.name=="Slc7a3")],nome_gene$Gene.stable.ID[which(nome_gene$Gene.name=="Pou5f1")])
geni_icb_ave=c(nome_gene$Gene.stable.ID[which(nome_gene$Gene.name=="Hhex")][1],nome_gene$Gene.stable.ID[which(nome_gene$Gene.name=="Cer1")],nome_gene$Gene.stable.ID[which(nome_gene$Gene.name=="Lefty1")])
geni_icb_ve=c(nome_gene$Gene.stable.ID[which(nome_gene$Gene.name=="Gata6")],nome_gene$Gene.stable.ID[which(nome_gene$Gene.name=="Apob")],nome_gene$Gene.stable.ID[which(nome_gene$Gene.name=="Afp")])
geni_icb=c(geni_icb_epi,geni_icb_exe,geni_icb_ve)
chosen.exprs_new_icb<-log10(data_tpm+1)[geni_icb,good_finale]
row_nuove<-nome_gene[row.names(chosen.exprs_new_icb),2]
geni_exe<-paste(c("Cdx2","Bmp4","Elf5"),"(ExE)",sep="")
geni_ep<-paste(c("Gng3","Slc7a3","Pou5f1"),"(Epi)",sep="")
geni_ve=paste(c("Gata6","Apob","Afp"),"(VE)",sep="")
row.names(chosen.exprs_new_icb)<-c(geni_ep,geni_exe,geni_ve)
numero_clu=rep(0,length(my.clusters))
numero_clu[cluster_finale==1]="Normal (Winner) Epiblast"
numero_clu[cluster_finale==2]="Extraembryonic Ectoderm"
numero_clu[cluster_finale==3]="Intermediate"
numero_clu[cluster_finale==4]="Loser Epiblast"
numero_clu[cluster_finale==5]="Visceral Endoderm"
#medium_clu=rep(0,length(medium_finale))
#medium_clu[medium_finale=="Cell competition OFF"]="CI-treated"
#medium_clu[medium_finale=="Cell competition ON"]="DMSO-treated"
haCol1 = HeatmapAnnotation(df = data.frame(Cluster = numero_clu
,Condition = medium_finale)
, col = list(Condition = c("Cell competition ON" = "yellow 1", "Cell competition OFF"="green 1"
)
,Cluster = c("Normal (Winner) Epiblast" = "#DD6400"
,"Extraembryonic Ectoderm" = "#800080"
,"Intermediate"= "#0000FF"
,"Loser Epiblast"="#006400"
,"Visceral Endoderm"="#56B4E9"
,"Losing "="#A9A9A9"
)
)
, show_legend = T
)
ht21 = Heatmap(chosen.exprs_new_icb# heat.vals_Mutant_Paired #
,cluster_rows = FALSE
, col= colorRamp2(c(0, round(max(chosen.exprs_new_icb))), c("white", "red"))
, name = "log 10 norm counts"
#, column_title = "Absolute values"
#, cluster_columns = as.dendrogram(my.tree)
,cluster_columns = F
#, cluster_rows =as.dendrogram(as.numeric(row.names(log_data_tpm)[1:10]) )
, top_annotation = haCol1
, row_names_gp = gpar(fontsize = 8
#,col = geneColors
)
, show_column_names = F
, show_row_names = T
)
## Warning: The input is a data frame, convert it to the matrix.
#setwd("/Users/gabriele.lubatti/Desktop/Phd/Cell_Competition/methdo_text_paper")
#pdf("heatmap_no_time.pdf",useDingbats = F)
draw(ht21
#,column_title = "Absolute values", column_title_side = "top"
)

#dev.off()
Extended Figure 2 C and D
row.names(classe)=good
cells_epi=row.names(classe)[classe$cluster==1|classe$cluster==3|classe$cluster==4]
Cluster_col=classe$cluster
library(ggplot2)
set.seed(1000)
dm=DiffusionMap(dissimilarity)
## Warning in DiffusionMap(dissimilarity): You have 723 genes. Consider passing
## e.g. n_pcs = 50 to speed up computation.
DC1=dm$DC1
DC2=dm$DC2
Cluster=as.factor(Cluster_col)
tre=data.frame(DC1,DC2,Cluster)
dpt <- DPT(dm)
#str(dpt)
#tips(dpt)
plot(dpt,col_by = 'DPT158')

cells_epi=row.names(classe)[classe$cluster==1|classe$cluster==3|classe$cluster==4]
set.seed(100)
dm=DiffusionMap(dissimilarity[cells_epi,cells_epi])
## Warning in DiffusionMap(dissimilarity[cells_epi, cells_epi]): You have 512
## genes. Consider passing e.g. n_pcs = 50 to speed up computation.
DC1=dm$DC1
DC2=dm$DC2
Cluster_col=classe[cells_epi,]$cluster
Cluster=as.factor(Cluster_col)
classe_epi=classe[cells_epi,]
Condition=as.factor(classe_epi$condizione)
tre=data.frame(DC1,DC2,Cluster)
diffu<-ggplot(tre, aes(x=DC1 , y=DC2 , color= Cluster)) + geom_point()
diffu + scale_color_manual(values=c("#DD6400", "#0000FF","#006400"),labels=c("Normal (winner) Epiblast","Intermediate", "Loser Epiblast"))+labs(x = "DC1",y="DC2") + ggtitle("Diffusion Map")

diffu<-ggplot(tre, aes(x=DC1 , y=DC2 , color= Condition)) + geom_point()+
scale_color_manual(values=c("green 1", "yellow 1"))
diffu

dpt <- DPT(dm)
#str(dpt)
#tips(dpt)
plot(dpt, col_by = 'DPT307')

time_all_epi=dpt$DPT307
Figure 2 E
cells_epi_fmk=row.names(classe)[(classe$cluster==1 & classe$condizione=="Cell competition OFF")|(classe$cluster==3 & classe$condizione=="Cell competition OFF")|(classe$cluster==4 & classe$condizione=="Cell competition OFF")]
cells_epi_dmso=row.names(classe)[(classe$cluster==1 & classe$condizione=="Cell competition ON")|(classe$cluster==3 & classe$condizione=="Cell competition ON")|(classe$cluster==4 & classe$condizione=="Cell competition ON")]
Cluster_clu_fmk=classe[cells_epi_fmk,]$cluster
Cluster_clu_dmso=classe[cells_epi_dmso,]$cluster
Condition_clu_fmk=classe[cells_epi_fmk,]$condizione
Condition_clu_dmso=classe[cells_epi_dmso,]$condizione
set.seed(1000)
dm=DiffusionMap(dissimilarity[cells_epi_fmk,cells_epi_fmk])
dmp=dissimilarity[cells_epi_dmso,cells_epi_fmk]
pred=dm_predict(dm,dmp)
prima=c(dm$DC1,pred[,1])
seconda=c(dm$DC2,pred[,2])
cluster=c(Cluster_clu_fmk,Cluster_clu_dmso)
Condition=c(Condition_clu_fmk,Condition_clu_dmso)
dati=data.frame(prima,seconda,cluster)
diffu<-ggplot(dati, aes(x=prima , y=seconda , color= as.factor(cluster))) + geom_point()
diffu + scale_color_manual(values=c("#DD6400", "#0000FF","#006400"),labels=c("Winning Epiblast","Losing Epiblast 2", "Losing Epiblast 1"))+labs(x = "DC1",y="DC2") + ggtitle("Diffusion Map")+labs(col="Cluster")

diffu<-ggplot(dati, aes(x=prima , y=seconda , color= as.factor(Condition))) + geom_point()
diffu<-ggplot(dati, aes(x=prima, y=seconda, color= Condition)) + geom_point()
#pdf("figure_paper_pdf/diffusion_map_epi_fmk_condition.pdf",width=8,height=5)
diffu + scale_color_manual(values=c("green 1", "yellow 1"))+labs(x = "DC1",y="DC2") + ggtitle("Diffusion Map")+labs(col="Condition")

dpt <- DPT(dm)
#save(dpt,file="dpt.Rda")
#str(dpt)
#tips(dpt)
#plot(dpt)
#plot(dpt,col=Cluster)
#pdf("figure_paper_pdf/diffusion_map_all_epi_pseudo_time.pdf",width=8,height=5)
plot(dpt, col_by = 'DPT119')

time_only_fmk=dpt$DPT119
projection_dist_ten=function (dm, new_dcs = NULL, ..., new_data, verbose = FALSE)
{
if (is.null(new_dcs))
new_dcs <- dm_predict(dm, new_data, ..., verbose = verbose)
else if (!missing(new_data))
stop("only pass one of new_dcs and new_data")
evs <- eigenvectors(dm)
nns <- find_knn(evs, 10, query = as.matrix(new_dcs))
nn_dist <- nns$index[, 1:10]
nn_dist
}
dist=projection_dist_ten(dm,new_dcs=pred)
k10=matrix(0,ncol=10,nrow = length(cells_epi_dmso))
for( i in 1:length(cells_epi_dmso)){
k10[i,]=time_only_fmk[dist[i,]]
#print(i)
}
time_projec=apply(k10,1,mean)
dfForPlot <- data.frame(clusters = as.character(Cluster_clu_fmk)
,exp=time_only_fmk
,Cluster = as.factor(Cluster_clu_fmk))
ggplot(data =dfForPlot
,aes(x = clusters
,y = exp
,col = Cluster )) +
geom_boxplot(width=0.9,alpha=0.8) +
scale_color_manual(breaks = c("1","3","4")
,labels = c("Normal (winner) Epiblast"
,"Intermediate"
,"Loser Epiblast "
)
,values=c("#DD6400","#0000FF","#006400" )) +
ylab("Losing Score")+xlab(NULL) +
geom_sina()+
ggtitle("Losing score FMK cells (274)")

dfForPlot <- data.frame(clusters = as.character(Cluster_clu_dmso)
,exp=time_projec
,Cluster = as.factor(Cluster_clu_dmso))
ggplot(data =dfForPlot
,aes(x = clusters
,y = exp
,col = Cluster )) +
geom_boxplot(width=0.9,alpha=0.8) +
scale_color_manual(breaks = c("1","3","4")
,labels = c("Normal (winner) Epiblast"
,"Intermediate"
,"Loser Epiblast "
)
,values=c("#DD6400","#0000FF","#006400" )) +
ylab("Losing Score")+xlab(NULL) +
geom_sina()+ylim(c(0,max(time_only_fmk)))+
ggtitle("Losing score DMSO cells (238)")

Figure 2F
primo=nome_gene$Gene.stable.ID[which(nome_gene$Gene.name=="Klf4")]
secondo=nome_gene$Gene.stable.ID[which(nome_gene$Gene.name=="Klf5")]
terzo=nome_gene$Gene.stable.ID[which(nome_gene$Gene.name=="Fgf5")]
quarto=nome_gene$Gene.stable.ID[which(nome_gene$Gene.name=="Zfp42")]
quinto=nome_gene$Gene.stable.ID[which(nome_gene$Gene.name=="Tdgf1")]
sesto=nome_gene$Gene.stable.ID[which(nome_gene$Gene.name=="Mesp1")]
settimo=nome_gene$Gene.stable.ID[which(nome_gene$Gene.name=="T")]
ottavo=nome_gene$Gene.stable.ID[which(nome_gene$Gene.name=="Neurod1")]
nono=nome_gene$Gene.stable.ID[which(nome_gene$Gene.name=="Sox1")]
decimo=nome_gene$Gene.stable.ID[which(nome_gene$Gene.name=="Sox17")]
quindici=nome_gene$Gene.stable.ID[which(nome_gene$Gene.name=="Gata6")]
undici=nome_gene$Gene.stable.ID[which(nome_gene$Gene.name=="Sox2")]
dodici=nome_gene$Gene.stable.ID[which(nome_gene$Gene.name=="Pou3f1")]
tredici=nome_gene$Gene.stable.ID[which(nome_gene$Gene.name=="Tcf7l1")]
quattordici=nome_gene$Gene.stable.ID[which(nome_gene$Gene.name=="Pou5f1")]
geni=c(primo,secondo,undici,dodici,tredici,quattordici,quarto,terzo,quinto,sesto,settimo,ottavo,nono,decimo,quindici)
chosen.exprs_new_icb<-log10(data_tpm+1)[geni,cells_epi_fmk]
#colnames(data_tpm)[order(time_all_tutto)]
new_versione=t(chosen.exprs_new_icb)[order(time_fmk_epi),]
chosen.exprs_new_icb=t(new_versione)
row_nuove<-nome_gene[row.names(chosen.exprs_new_icb),2]
#row.names(chosen.exprs_new_icb)<-row_nuove
geni_naive<-paste(c("Klf4*","Klf5*"),"(Naive)",sep="")
geni_epi<-paste(c("Fgf5*","Tdgf1*"),"(Epiblast)",sep="")
#geni_ave<-paste(c("Hhex","Cer1","Lefty1"),"(AVE)",sep="")
geni_meso=paste(c("Mesp1","T"),"(Mesoderm)",sep="")
geni_neuro=paste(c("Neurod1*","Sox1"),"(Neuroderm)",sep="")
geni_endo=paste(c("Sox17*","Gata6"),"(Endoderm)",sep="")
geni_new=paste(c("Sox2","Pou3f1","Tcf7l1*","Pou5f1*","Rex1"),"(Naive)",sep="")
row.names(chosen.exprs_new_icb)<-c(geni_naive,geni_new,geni_epi,geni_meso,geni_neuro,geni_endo)
classe_fmk=meta_dati_course[cells_epi_fmk,]
classe_sort=classe_fmk[order(time_fmk_epi),]
new_clusters=classe_sort$cluster
Condition_new=as.vector(classe_sort$condition)
numero_clu=rep(0,length(new_clusters))
numero_clu[new_clusters==1]="Normal (Winner) Epiblast"
numero_clu[new_clusters==2]="Extraembryonic Ectoderm"
numero_clu[new_clusters==3]="Intermediate"
numero_clu[new_clusters==4]="Loser Epiblast"
numero_clu[new_clusters==5]="Visceral Endoderm"
numero_clu[new_clusters==6]="Losing "
bo=data.frame(Condition = Condition_new,Cluster = numero_clu,
Pseudotime=sort(time_fmk_epi))
haCol1 = HeatmapAnnotation(df = data.frame(Condition = Condition_new
,Cluster = numero_clu,
Pseudotime=sort(time_fmk_epi))
, col = list(Condition = c("Cell competition ON" = "yellow 1", "Cell competition OFF"="green 1"
)
,Cluster = c("Normal (Winner) Epiblast" = "#DD6400"
,"Extraembryonic Ectoderm" = "#800080"
,"Intermediate"= "#0000FF"
,"Loser Epiblast"="#006400"
,"Visceral Endoderm"="#56B4E9"
,"Losing "="#A9A9A9"
)
,Pseudotime=colorRamp2(c(0, round(max(time_fmk_epi))), c("white", "blue")))
, show_legend = T
)
#pdf("figure_paper_pdf/heatmap_pseudo_time_naive.pdf",width=8,height=5)
ht21 = Heatmap(chosen.exprs_new_icb# heat.vals_Mutant_Paired #
,cluster_rows = FALSE
, col= colorRamp2(c(0, round(max(chosen.exprs_new_icb))), c("white", "red"))
, name = "log 10 norm counts"
#, column_title = "Absolute values"
, cluster_columns = F
,column_title = " Naive and primitive streak markers"
#, cluster_rows =as.dendrogram(as.numeric(row.names(log_data_tpm)[1:10]) )
, top_annotation = haCol1
, row_names_gp = gpar(fontsize = 8
#,col = geneColors
)
, show_column_names = F
, show_row_names = T
)
draw(ht21
#,column_title = "Absolute values", column_title_side = "top"
)

Extended Figure 3 D
singolo_correlazione=function(a,b,c,d){
classe_prova=classe[grep(b,classe$Cell_type),]
classe_los=classe_prova[grep(a,classe_prova$cluster),]
exp_1_los=t((data_tpm))[rownames(classe_los),nome_gene$Gene.stable.ID[which(nome_gene$Gene.name==c)]]
exp_1_los=as.vector(exp_1_los)
exp_1_los_1=t((data_tpm))[rownames(classe_los),nome_gene$Gene.stable.ID[which(nome_gene$Gene.name==d)]]
exp_1_los_1=as.vector(exp_1_los_1)
dati=data.frame(exp_1_los,exp_1_los_1)
p=ggplot(data=dati,aes(x=exp_1_los_1,y=exp_1_los))+geom_point(col="#006400")+xlab(d)+
ylab(c)+ggtitle(paste(c,"vs",d,"-spearman corr:",round(cor(exp_1_los_1, exp_1_los, method = "spearman"),3),sep=" "))
p
}
singolo_correlazione("4","FMK","Klf4","Sox17")

singolo_correlazione("4","FMK","Klf4","Neurod1")

singolo_correlazione("4","FMK","Neurod1","Sox17")

GAM Model for detecting DE genes along the pseudo time using only Cell Competition OFF epiblast cells
row.names(classe)=good
cellule_epi_fmk=row.names(classe)[(classe$cluster==1 & classe$condizione=="Cell competition OFF")|(classe$cluster==3 & classe$condizione=="Cell competition OFF")|(classe$cluster==4 & classe$condizione=="Cell competition OFF")]
cellule_epi=row.names(classe)[(classe$cluster==1)|(classe$cluster==3) |(classe$cluster==4 )]
data=data_tpm
data=data[,cellule_epi_fmk]
risultato=apply(data,1,function(x){sum(x>15)})
indice=which(as.vector(risultato)>10)
background_trajectory=row.names(data_tpm)[indice]
data=data[indice,]
data=log10(data+1)
sce<-SingleCellExperiment(assays = list(counts = as.matrix(data)))
sce$Pseudotime <- time_fmk_epi
Y=counts(sce)
#Y_1=Y[1:2,]
t=sce$Pseudotime
gam.res <- apply(Y, 1, function(z){
d <- data.frame(z=z, t=t)
tmp <- gam(z ~ lo(t), data=d)
p <- summary(tmp)[4][[1]][1,5] #p-value parametric ANOVA
f<-fitted(tmp)
#p2 <- summary(tmp)[3][[1]][2,3] #p-value non-parametric ANOVA
c(p,f)
#c(p1,p2) #If we want to return both the parametric and non-parametric p-values
})
data_all=data
#d <- data.frame(z=Y[1,], t=t)
#tmp <- gam(z ~ lo(t), data=d)
#p <- summary(tmp)[4][[1]][1,5] #p-value parametric ANOVA
#f<-fitted(tmp)
Keeping only the genes which vary in a significant way along the pseudo time
genes.table<-data.frame(genes.names=rownames(data))
genes.table$pvals<-gam.res[1,]
genes.table$FDR<-p.adjust(genes.table$pvals, method="fdr") #Adjust p-values
genes.table<-genes.table[order(genes.table$FDR),] #Order the genes according to the FDR
genes.table$genes.names <- as.character(genes.table$genes.names)
#Note that we are using a very strict threshold on FDR! Otherwise we obtain a very large number
#of genes
results.gam.tot<-genes.table[genes.table$FDR < 0.01,][c('genes.names','FDR')] #Filter the significant genes
#Save the matrix of fitted values
gam.fitted<-gam.res[-1,]
print(paste("We have"
,sum(genes.table$FDR < 0.01)
,"significant genes"))
## [1] "We have 5310 significant genes"
#results_gam=results.gam.tot['genes.names']
results_gam_nomi=results.gam.tot$genes.names
gam_fitted=gam.fitted
data=data[results_gam_nomi,]
data=t(data)
data_100_sort=data
Hierarchical cluster analysis for detecting genes who behave in a similar way along the trajectory
chosen.exprs1<-as.matrix(data_100_sort)
dissimilarity <- sqrt((1-cor((chosen.exprs1), method = "spearman"))/2)
#Umap by medium
my.dist <- as.dist(dissimilarity)
set.seed(100)
#object<-umap(d=as.matrix(my.dist),input="dist")
#save(object,file="object_novembre.Rda")
object<-umap(d=as.matrix(my.dist),input="dist")
umap_met<-as.data.frame(object$layout)
#umap_met
my.tree <- hclust(my.dist
, method="ward.D2" # met = "ward.D2" was adopted form the reference
)
# Cut the treee to identify clusters
#save(my.tree,file="my_tree_novembre.Rda")
my.clusters <- cutreeHybrid(my.tree
,distM=as.matrix(my.dist),deepSplit=0
,minClusterSize = 100)$label#5)$label
## ..cutHeight not given, setting it to 14.3 ===> 99% of the (truncated) height range in dendro.
## ..done.
#save(my.clusters,file="my_clusters_novembre.Rda")
#length(unique(my.clusters))
#classe<-data.frame(Cell_type = meta_dati_no$Sample
#,cluster = my.clusters)
#classe$condizione=Condition
#table(classe$cluster,classe$condizione)
Cluster=as.factor(my.clusters)
sp<-ggplot(as.data.frame(object$layout),aes(umap_met[,1],umap_met[,2],color=Cluster)) + geom_point()+xlab("Umap 1")+ylab("Umap 2")+ ggtitle("Cluster on top (FDR<0.01) varying genes along pseudotime")
sp

Identification of DE genes (up= increasing along the trajectory, down= decreasing along the trajectory)
classe_geni=data.frame(colnames(data_100_sort),my.clusters)
colnames(classe_geni)=c("geni","clusters")
row.names(classe_geni)=classe_geni$geni
geni_1=nome_gene[(classe_geni[classe_geni$clusters==1,]$geni),2]
geni_2=nome_gene[(classe_geni[classe_geni$clusters==2,]$geni),2]
geni_3=nome_gene[(classe_geni[classe_geni$clusters==3,]$geni),2]
geni_4=nome_gene[(classe_geni[classe_geni$clusters==4,]$geni),2]
id_1=as.vector(classe_geni$geni[classe_geni$clusters==1])
id_2=as.vector(classe_geni$geni[classe_geni$clusters==2])
id_3=as.vector(classe_geni$geni[classe_geni$clusters==3])
id_4=as.vector(classe_geni$geni[classe_geni$clusters==4])
id_down_loser=c(id_1,id_2,id_4)
id_up_loser=c(id_3)
row.names(results.gam.tot)=results.gam.tot$genes.names
id_1_full=results.gam.tot[id_1,]
id_1_full$Cluster=classe_geni[id_1,]$clusters
id_2_full=results.gam.tot[id_2,]
id_2_full$Cluster=classe_geni[id_2,]$clusters
id_3_full=results.gam.tot[id_3,]
id_3_full$Cluster=classe_geni[id_3,]$clusters
id_4_full=results.gam.tot[id_4,]
id_4_full$Cluster=classe_geni[id_4,]$clusters
id_down_full=results.gam.tot[id_down_loser,]
id_down_full$Cluster=classe_geni[id_down_loser,]$clusters
id_down_full_sort=id_down_full[order(id_down_full$FDR),]
id_down_full_sort$Cluster=classe_geni[row.names(id_down_full_sort),]$clusters
Extended Figure 3 A
row.names(classe)=good
cellule_epi_fmk=row.names(classe)[(classe$cluster==1 & classe$condizione=="Cell competition OFF")|(classe$cluster==3 & classe$condizione=="Cell competition OFF")|(classe$cluster==4 & classe$condizione=="Cell competition OFF")]
geni_up=results_gam_nomi[results_gam_nomi%in%id_up_loser]
geni_down=results_gam_nomi[results_gam_nomi%in%id_down_loser]
geni=c(geni_up,geni_down)
geni_colore=rep(0,length(geni))
geni_colore[1:length(geni_up)]="Increasing"
geni_colore[(length(geni_up)+1):length(geni)]="Decreasing"
DE_genes=geni_colore
chosen.exprs_new_icb<-log10(data_tpm+1)[geni,cellule_epi_fmk]
new_versione=t(chosen.exprs_new_icb)[order(time_fmk_epi),]
chosen.exprs_new_icb=t(new_versione)
row_nuove<-nome_gene[row.names(chosen.exprs_new_icb),2]
geni_nomi=nome_gene[geni,]$Gene.name
row.names(chosen.exprs_new_icb)<-geni_nomi
row.names(classe)=good
classe_fmk=classe[cellule_epi_fmk,]
classe_sort=classe_fmk[order(time_fmk_epi),]
new_clusters=classe_sort$cluster
Condition_new=classe_sort$condizione
numero_clu=rep(0,length(new_clusters))
numero_clu[new_clusters==1]="Normal (Winner) Epiblast"
numero_clu[new_clusters==2]="Extraembryonic Ectoderm"
numero_clu[new_clusters==3]="Intermediate"
numero_clu[new_clusters==4]="Loser Epiblast"
numero_clu[new_clusters==5]="Visceral Endoderm"
numero_clu[new_clusters==6]="Losing "
haCol1 = HeatmapAnnotation(df = data.frame(
Cluster = numero_clu,
Condition = Condition_new,
Pseudotime=sort(time_fmk_epi))
, col = list(Condition = c("Cell competition ON" = "yellow 1", "Cell competition OFF"="green 1"
)
,Cluster = c("Normal (Winner) Epiblast" = "#DD6400"
,"Extraembryonic Ectoderm" = "#800080"
,"Intermediate"= "#0000FF"
,"Loser Epiblast"="#006400"
,"Visceral Endoderm"="#56B4E9"
,"Losing "="#A9A9A9"
)
,Pseudotime=colorRamp2(c(0, round(max(time_fmk_epi))), c("white", "blue")))
, show_legend = T,which = c("column")
)
left_annotation=HeatmapAnnotation(df = data.frame(
DE_genes=DE_genes)
, col = list(DE_genes = c("Increasing" = "black", "Decreasing"="grey")
)
, show_legend = T,which = c("row")
)
ht21 = Heatmap(chosen.exprs_new_icb# heat.vals_Mutant_Paired #
,cluster_rows = FALSE
, col= colorRamp2(c(0,round(max(chosen.exprs_new_icb))), c("blue","red"))
, name = "log 10 norm counts"
#, column_title = "Absolute values"
, cluster_columns = F
,column_title = "Diff expressed genes"
#, cluster_rows =as.dendrogram(as.numeric(row.names(log_data_tpm)[1:10]) )
, top_annotation = haCol1
, row_names_gp = gpar(fontsize = 8
#,col = geneColors
)
, show_column_names = F
, show_row_names = F
,left_annotation = left_annotation )
## `use_raster` is automatically set to TRUE for a matrix with more than
## 2000 rows. You can control `use_raster` arugment by explicitly setting
## TRUE/FALSE to it. Set `ht_opt$message = FALSE` to turn off this
## message.
draw(ht21
#,column_title = "Absolute values", column_title_side = "top"
)

Extended Figure 3 B and C
up_loser_1_4=data.frame(id_up_loser)
down_loser_1_4=data.frame(id_down_loser)
row.names(up_loser_1_4)=id_up_loser
row.names(down_loser_1_4)=id_down_loser
target_up=data_p53$Ensembl.ID[grep("up",data_p53$Literature..Table.S1.)]
target_down=data_p53$Ensembl.ID[grep("down",data_p53$Literature..Table.S1.)]
#write.csv2(target_up,file="target_up_p53.csv")
#write.csv2(target_down,file="target_down_p53.csv")
#ho generato i seguenti due file con g profiler
#up_p53_mouse=read.csv("up_p53_topo.csv",sep=",",header=T)
#down_p53_mouse=read.csv("down_p53_topo.csv",sep=",",header=T)
#test up
test=up_p53_mouse$initial_ensg
test= !(duplicated(test) | duplicated(test, fromLast = TRUE)) # correct result
up_p53=up_p53_mouse$ortholog_ensg[test]
up_p53=up_p53[up_p53%in%row.names(down_loser_1_4)|up_p53%in%row.names(up_loser_1_4)]
sum(up_p53%in%row.names(up_loser_1_4))/(length(up_p53))
## [1] 0.5625
sum(up_p53%in%row.names(down_loser_1_4))/(length(up_p53))
## [1] 0.4375
#test down
test=down_p53_mouse$initial_ensg
test= !(duplicated(test) | duplicated(test, fromLast = TRUE)) # correct result
down_p53=down_p53_mouse$ortholog_ensg[test]
down_p53=down_p53[down_p53%in%row.names(down_loser_1_4)|down_p53%in%row.names(up_loser_1_4)]
sum(down_p53%in%row.names(up_loser_1_4))/(length(down_p53))
## [1] 0.195122
sum(down_p53%in%row.names(down_loser_1_4))/(length(down_p53))
## [1] 0.804878
tot_p53=c(up_p53,down_p53)
matrice=matrix(c(sum(up_p53%in%row.names(up_loser_1_4)),length(up_p53)-sum(up_p53%in%row.names(up_loser_1_4)),sum(down_p53%in%row.names(up_loser_1_4)),sum(down_p53%in%row.names(down_loser_1_4))),nrow=2,ncol=2,byrow=T)
#Fisher test
fisher.test(matrice)
##
## Fisher's Exact Test for Count Data
##
## data: matrice
## p-value = 0.0001076
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 2.046214 14.816715
## sample estimates:
## odds ratio
## 5.229645
matrice_1=matrix(c(sum(down_p53%in%row.names(down_loser_1_4)),length(down_p53)-sum(down_p53%in%row.names(down_loser_1_4)),sum(up_p53%in%row.names(down_loser_1_4)),sum(up_p53%in%row.names(up_loser_1_4))),nrow=2,ncol=2,byrow=T)
fisher.test(matrice_1)
##
## Fisher's Exact Test for Count Data
##
## data: matrice_1
## p-value = 0.0001076
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 2.046214 14.816715
## sample estimates:
## odds ratio
## 5.229645
a=rep("Up genes (45)",1)
b=rep("Down genes (35)",1)
aa=sum(up_p53%in%row.names(up_loser_1_4))
bb=sum(up_p53%in%row.names(down_loser_1_4))
value=c(aa,bb)
group=c(a,b)
df=data.frame(group,value)
library(ggplot2)
bp<- ggplot(df, aes(x="", y=value, fill=group))+
geom_bar(width = 1, stat = "identity")
#bp
pie <- bp + coord_polar("y", start=0)
#pie
library(scales)
##
## Attaching package: 'scales'
## The following object is masked from 'package:viridis':
##
## viridis_pal
blank_theme <- theme_minimal()+
theme(
axis.title.x = element_blank(),
axis.title.y = element_blank(),
panel.border = element_blank(),
panel.grid=element_blank(),
axis.ticks = element_blank(),
plot.title=element_text(size=14, face="bold")
)
pie + blank_theme +
theme(axis.text.x=element_blank()) +
geom_text(aes(y=c(6,60),
label = percent(value/80)), size=5)+ggtitle("Targets activated by P53")

a=rep("Up genes (8)",1)
b=rep("Down genes (33)",1)
aa=sum(down_p53%in%row.names(up_loser_1_4))
bb=sum(down_p53%in%row.names(down_loser_1_4))
value=c(aa,bb)
group=c(a,b)
df=data.frame(group,value)
library(ggplot2)
bp<- ggplot(df, aes(x="", y=value, fill=group))+
geom_bar(width = 1, stat = "identity")
#bp
pie <- bp + coord_polar("y", start=0)
#pie
library(scales)
blank_theme <- theme_minimal()+
theme(
axis.title.x = element_blank(),
axis.title.y = element_blank(),
panel.border = element_blank(),
panel.grid=element_blank(),
axis.ticks = element_blank(),
plot.title=element_text(size=14, face="bold")
)
pie + blank_theme +
theme(axis.text.x=element_blank()) +
geom_text(aes(y=c(4,30),
label = percent(value/41)), size=5)+ggtitle("Targets repressed by P53")

gene_expression <- function(a,b){
#row.names(clas_new_new)=good
geni_up_loser=row.names(classe_geni)[classe_geni$clusters==3]
geni_down_loser=row.names(classe_geni)[classe_geni$clusters==1|classe_geni$clusters==2|classe_geni$clusters==4]
classe_new=classe[classe$cluster==1|classe$cluster==3|classe$cluster==4,]
clas_new_new=classe_new
clas_new_new=clas_new_new[grep("FMK",classe_new$Cell_type),]
classe=clas_new_new
exp_real=t(log10(data_tpm+1))[rownames(classe),nome_gene$Gene.stable.ID[which(nome_gene$Gene.name==a)]]
#exp_fitted=gam.fitted[rownames(classe),which(colnames(gam.fitted)==nome_gene$Gene.stable.ID[which(nome_gene$Gene.name==a)])]
#exp=c(exp_real)
pseudo_time_real=time_fmk_epi
#pseudo_time_fit=time_all
#type=c(rep("real expression",length(exp_real)),rep("fitted value",length(exp_real)))
#time_all_all=c(pseudo_time_real,pseudo_time_fit)
dfForPlot <- data.frame(pseudo_time = pseudo_time_real
,genes_expression=exp_real
)
if(sum(nome_gene$Gene.stable.ID[which(nome_gene$Gene.name==a)]%in%geni_up_loser)>0){
titolo=paste(a,"up_loser",b,sep=" ")
}
if(sum(nome_gene$Gene.stable.ID[which(nome_gene$Gene.name==a)]%in%geni_down_loser)>0){
titolo=paste(a,"down_loser",b,sep=" ")
}
if(sum(nome_gene$Gene.stable.ID[which(nome_gene$Gene.name==a)]%in%colnames(gam.fitted))==0){titolo=paste(a,"no in the background",b,sep=" ")}
if(sum((nome_gene$Gene.stable.ID[which(nome_gene$Gene.name==a)]%in%colnames(gam.fitted)>0)&(nome_gene$Gene.stable.ID[which(nome_gene$Gene.name==a)]%in%row.names(classe_geni)==0))){titolo=paste(a,"no up/down",b,sep=" ")}
if(sum(nome_gene$Gene.stable.ID[which(nome_gene$Gene.name==a)]%in%colnames(gam.fitted))==0){titolo=paste(a,"no in the background",b,sep=" ")}
p <- ggplot(data =dfForPlot
,aes(x = pseudo_time
,y = genes_expression
)) +
geom_point()+xlab("pseudotime")+ylab("log 10 norm counts")+
ggtitle(titolo)+geom_smooth(method = "loess", se = F,col="red")
p}
Figure 3 C and Extended Figure 4 D
p_12=gene_expression("mt-Atp6"," ")
p_13=gene_expression("mt-Nd3"," ")
p_14=gene_expression("Opa1"," ")
p_15=gene_expression("Samm50"," ")
#pdf("figure_paper_pdf/pseudo_time_elctron.pdf",width=8,height=5)
grid.arrange(p_12,p_13,top=textGrob("Electron Transport chain",gp=gpar(fontsize=15)),ncol=1,nrow=2)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

grid.arrange(p_14,p_15,top=textGrob("Mitochondria Membrane",gp=gpar(fontsize=15)),ncol=1,nrow=2)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

p_12=gene_expression("Atf4"," ")
p_13=gene_expression("Ddit3"," ")
p_14=gene_expression("Foxo3"," ")
p_15=gene_expression("Nfe2l2"," ")
p_16=gene_expression("Sesn2"," ")
grid.arrange(p_14,p_15,ncol=1,nrow=2)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

grid.arrange(p_12,p_13,ncol=1,nrow=2)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

grid.arrange(p_16,ncol=1,nrow=1)
## `geom_smooth()` using formula 'y ~ x'

gene_expression <- function(a,b){
#row.names(clas_new_new)=good
geni_up_loser=row.names(classe_geni)[classe_geni$clusters==3]
geni_down_loser=row.names(classe_geni)[classe_geni$clusters==1|classe_geni$clusters==2|classe_geni$clusters==4]
classe_new=classe[classe$cluster==1|classe$cluster==3|classe$cluster==4,]
clas_new_new=classe_new
clas_new_new=clas_new_new[grep("FMK",classe_new$Cell_type),]
classe=clas_new_new
cluster_plot=classe[,2]
cluster_plot[cluster_plot==1]="Winner Epiblast"
cluster_plot[cluster_plot==3]="Intermediate"
cluster_plot[cluster_plot==4]="Loser Epiblast"
cluster_plot=factor(cluster_plot,levels=c("Winner Epiblast","Intermediate","Loser Epiblast"))
exp_real=t(log10(data_tpm+1))[rownames(classe),nome_gene$Gene.stable.ID[which(nome_gene$Gene.name==a)]]
#exp_fitted=gam.fitted[rownames(classe),which(colnames(gam.fitted)==nome_gene$Gene.stable.ID[which(nome_gene$Gene.name==a)])]
#exp=c(exp_real)
pseudo_time_real=time_fmk_epi
#pseudo_time_fit=time_all
#type=c(rep("real expression",length(exp_real)),rep("fitted value",length(exp_real)))
#time_all_all=c(pseudo_time_real,pseudo_time_fit)
dfForPlot <- data.frame(pseudo_time = pseudo_time_real
,genes_expression=exp_real
,cluster=cluster_plot
)
if(sum(nome_gene$Gene.stable.ID[which(nome_gene$Gene.name==a)]%in%geni_up_loser)>0){
titolo=paste(a,"up_loser",b,sep=" ")
}
if(sum(nome_gene$Gene.stable.ID[which(nome_gene$Gene.name==a)]%in%geni_down_loser)>0){
titolo=paste(a,"down_loser",b,sep=" ")
}
if(sum(nome_gene$Gene.stable.ID[which(nome_gene$Gene.name==a)]%in%colnames(gam.fitted))==0){titolo=paste(a,"no in the background",b,sep=" ")}
if(sum((nome_gene$Gene.stable.ID[which(nome_gene$Gene.name==a)]%in%colnames(gam.fitted)>0)&(nome_gene$Gene.stable.ID[which(nome_gene$Gene.name==a)]%in%row.names(classe_geni)==0))){titolo=paste(a,"no up/down",b,sep=" ")}
if(sum(nome_gene$Gene.stable.ID[which(nome_gene$Gene.name==a)]%in%colnames(gam.fitted))==0){titolo=paste(a,"no in the background",b,sep=" ")}
p <- ggplot(data =dfForPlot
,aes(x = pseudo_time
,y = genes_expression
,color=(cluster_plot)
)) +
geom_point()+xlab("pseudotime")+ylab("log 10 norm counts")+
ggtitle(titolo)+geom_smooth(method = "loess", se = F,col="red") +scale_color_manual(values=c("#DD6400", "#0000FF", "#006400"))+labs(col="Cluster")
p <- ggplot(data =dfForPlot
,aes(x = cluster_plot
,y = genes_expression
,colour=(cluster_plot)
)) +
geom_boxplot()+xlab("Cluster")+ylab("log 10 norm counts")+
ggtitle(titolo) +scale_colour_manual(values=c("#DD6400", "#0000FF", "#006400"))+labs(colour="Cluster")
list(p,q)}
#("Normal (Winner) Epiblast" = "#DD6400"
# ,"Extraembryonic Ectoderm" = "#800080"
#,"Intermediate"= "#0000FF"
# ,"Loser Epiblast"="#006400"
# ,"Visceral Endoderm"="#56B4E9"
# ,"Losing "="#A9A9A9"
#gene_expression("Tfb1m"," ")
#setwd("/Users/gabriele.lubatti/Desktop")
#geni_upr=read.csv2("ana_upr.csv",header=F)
geni_upr=as.vector(geni_upr$V1)
geni_down=nome_gene[id_down_loser,]$Gene.name
geni_up=nome_gene[id_up_loser,]$Gene.name
geni_upr[which(!(geni_upr%in%nome_gene$Gene.name))]="Foxo3"
Extended Figure 4 B and C
geni_up=nome_gene[row.names(id_3_full),]$Gene.name
id_3_full$nomi=geni_up
geni_up_com=geni_up[which(geni_up%in%geni_upr)]
index_up_com=which(geni_up%in%geni_upr)
fdr_up_com=id_3_full[which(geni_up%in%geni_upr),]$FDR
up_com_1=data.frame(geni_up_com,fdr_up_com,index_up_com)
up_com_1
## geni_up_com fdr_up_com index_up_com
## 1 Ddit3 4.746395e-39 2
## 2 Atf3 6.102567e-27 22
## 3 Atf4 2.201811e-23 31
## 4 Foxo3 2.738259e-22 37
## 5 Ppp1r15a 8.494053e-18 68
## 6 Eif2ak3 7.200497e-13 150
## 7 Nfe2l2 1.563468e-10 207
## 8 Gdf15 5.512919e-08 331
geni_down=nome_gene[row.names(id_down_full_sort),]$Gene.name
id_down_full_sort=id_down_full_sort[,-1]
id_down_full_sort$nomi=geni_down
geni_down_com=geni_down[which(geni_down%in%geni_upr)]
index_down_com=which(geni_down%in%geni_upr)
fdr_down_com=id_down_full_sort[which(geni_down%in%geni_upr),]$FDR
down_com_1=data.frame(geni_down_com,fdr_down_com,index_down_com)
down_com_1
## geni_down_com fdr_down_com index_down_com
## 1 Mthfd1l 2.504068e-35 147
## 2 Hspe1 8.853980e-34 164
## 3 Cat 2.428744e-30 219
## 4 Hspd1 6.863457e-13 1262
## 5 Sod2 1.255766e-10 1552
## 6 Hsph1 4.447044e-10 1654
## 7 Lonp1 1.061746e-06 2346
## 8 Eif2a 1.469514e-06 2381
## 9 Mthfd2 1.297814e-05 2691
## 10 Hspa4 2.819542e-05 2789
## 11 Cth 2.533223e-03 3677
## 12 Nrf1 2.843621e-03 3695
back=background_trajectory
#versione corretta con geni id
id_geni_upr=nome_gene$Gene.stable.ID[nome_gene$Gene.name%in%geni_upr]
#arricchimento per up
matrice=matrix(c(sum(id_up_loser%in%id_geni_upr),sum(!id_up_loser%in%id_geni_upr),sum(back[!(back%in%id_up_loser)]%in%id_geni_upr),sum(!(back[!(back%in%id_up_loser)]%in%id_geni_upr))),nrow=2,ncol=2,byrow = T)
fisher.test(matrice)
##
## Fisher's Exact Test for Count Data
##
## data: matrice
## p-value = 0.01271
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 1.148861 7.414376
## sample estimates:
## odds ratio
## 3.06029
#arricchimento per down
matrice=matrix(c(sum(id_down_loser%in%id_geni_upr),sum(!(id_down_loser%in%id_geni_upr)),sum(back[!(back%in%id_down_loser)]%in%id_geni_upr),sum(!(back[!(back%in%id_down_loser)]%in%id_geni_upr))),nrow=2,ncol=2,byrow = T)
fisher.test(matrice)
##
## Fisher's Exact Test for Count Data
##
## data: matrice
## p-value = 0.6902
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.524691 2.899567
## sample estimates:
## odds ratio
## 1.243201
setwd(base_path)
raw_dati=read.table("GeneCounts.csv",sep=",",header=T)
row.names(raw_dati)=raw_dati$X
raw_dati=raw_dati[,-1]
colnames(raw_dati)[9]="P2A4HB_GFP_CO_1"
colnames(raw_dati)[10]="P2A4HB_GFP_CO_2"
#INIZIO ARRAY EXPRESS
setwd(base_path)
meta_data=read.table("metadata_ana_new.csv",sep=";",header=F)
nuovo=as.vector(meta_data$V5)
nuovo[meta_data$V5=="Separate"]="Sep"
nuovo[meta_data$V5=="Co-culture"]="CO"
meta_data$V5=nuovo
nomi_colonne=paste(meta_data$V4,meta_data$V5,meta_data$V6,sep="_")
nomi_colonne[which(!nomi_colonne%in%colnames(raw_dati))]=c("P2A4HB_GFP_Sep_1", "P2A4HB_GFP_CO_1", "P2A4HB_GFP_Sep_2", "P2A4HB_GFP_CO_2","P2A4HB_GFP_Sep_3", "P2A4HB_GFP_CO_3", "P2A4HB_GFP_Sep_4", "P2A4HB_GFP_CO_4" )
nuovi_dati=raw_dati[,nomi_colonne]
colnames(nuovi_dati)=as.vector(meta_data$V3)
De_edger=function(Treat,Data){
Treat=relevel(Treat, ref="Win_Epi")
y=DGEList(counts=Data, group=Treat)
keep <- filterByExpr(y)
tutto=keep
table(keep)
y=y[keep, , keep.lib.sizes=FALSE]
y=calcNormFactors(y)
#y$samples
design <- model.matrix(~Treat)
rownames(design) <- colnames(y)
#design
y=estimateDisp(y, design, robust=TRUE)
#y$common.dispersion
fit=glmQLFit(y, design, robust=TRUE)
#qlf=glmQLFTest(fit, coef=1:5)
#qlf=glmLRT(fit)
#topTags(qlf)
#FDR= p.adjust(qlf$table$PValue, method="BH")
#qlf$table<-cbind(qlf$table,FDR)
#sum(FDR < 0.05)
qlf=glmQLFTest(fit)
FDR= p.adjust(qlf$table$PValue, method="BH")
qlf$table<-cbind(qlf$table,FDR)
sum(FDR < 0.05)
topTags(qlf)
summary(decideTests(qlf))
#plotMD(qlf)
#abline(h=c(-1,1), col="blue")
tr_new<-qlf$table[row.names(qlf$table)[qlf$table$FDR<0.001 & !is.na(qlf$table$FDR)],]
#return(list(tr_new))
#}
tr_export_down<-tr_new[row.names(tr_new)[tr_new$logFC<0],]
tr_down<-tr_export_down[order(tr_export_down$FDR),]
tr_down_1_4_new<-tr_down
tr_export_up<-tr_new[row.names(tr_new)[tr_new$logFC>0],]
tr_up<-tr_export_up[order(tr_export_up$FDR),]
tr_up_1_4_new<-tr_up
name_up=nome_gene[row.names(tr_up_1_4_new),2]
tr_up_1_4_new$name=name_up
name_down=nome_gene[row.names(tr_down_1_4_new),2]
tr_down_1_4_new$name=name_down
return(list(tr_up_1_4_new,tr_down_1_4_new,tutto))
}
#prima funzione
colnames(raw_dati)[9]="P2A4HB_GFP_CO_1"
colnames(raw_dati)[10]="P2A4HB_GFP_CO_2"
winner=colnames(raw_dati)[grep("P2A4HB_GFP_CO",colnames(raw_dati))]
loser=colnames(raw_dati)[grep("P1A3BG_CO",colnames(raw_dati))]
cel=c(winner,loser)
data_de=raw_dati[,cel]
tre_1=(rep("Win_Epi",length(winner)))
tre_2=(rep("Los_Epi",length(loser)))
Treat=c(tre_1,tre_2)
Treat=as.factor(Treat)
#Time=as.vector((classe[cel,]$Cell_type))
#Time=as.factor(Time)
Output_De=De_edger(Treat,data_de)
up_1_4=Output_De[[1]]
down_1_4=Output_De[[2]]
tutto=Output_De[[3]]
background_all=row.names(raw_dati)[tutto]
#id_down_loser
#id_up_loser
back_all=intersect(background_trajectory,background_all)
Extended Figure 10 B
library(UpSetR)
up_1_3_n=id_up_loser
down_1_3_n=id_down_loser
up_1_4_n=row.names(up_1_4)
down_1_4_n=row.names(down_1_4)
listInput <- list(up_vivo = up_1_3_n, down_vivo =down_1_3_n, up_vitro = up_1_4_n,down_vitro=down_1_4_n)
upset(fromList(listInput), order.by = "freq")

down_1_4=down_1_4[intersect(row.names(down_1_4),back_all),]
up_1_4=up_1_4[intersect(row.names(up_1_4),back_all),]
id_up_loser=intersect(id_up_loser,back_all)
id_down_loser=intersect(id_down_loser,back_all)
#arrichimento per down
matrice=matrix(c(sum(row.names(down_1_4)%in%id_down_loser),sum(row.names(down_1_4)%in%back_all)-sum(row.names(down_1_4)%in%id_down_loser),sum(id_down_loser%in%back_all)-sum(row.names(down_1_4)%in%id_down_loser),sum(!(back_all%in%row.names(down_1_4))&!(back_all%in%id_down_loser))),nrow=2,ncol=2,byrow = T)
#Fisher test
fisher.test(matrice)
##
## Fisher's Exact Test for Count Data
##
## data: matrice
## p-value = 1.706e-12
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 1.523194 2.128908
## sample estimates:
## odds ratio
## 1.799954
#arrichimento per up
matrice=matrix(c(sum(row.names(up_1_4)%in%id_up_loser),sum(row.names(up_1_4)%in%back_all)-sum(row.names(up_1_4)%in%id_up_loser),sum(id_up_loser%in%back_all)-sum(row.names(up_1_4)%in%id_up_loser),sum(!(back_all%in%row.names(up_1_4))&!(back_all%in%id_up_loser))),nrow=2,ncol=2,byrow = T)
#Fisher test
fisher.test(matrice)
##
## Fisher's Exact Test for Count Data
##
## data: matrice
## p-value = 0.3314
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.8780489 1.4101126
## sample estimates:
## odds ratio
## 1.117812
#arrichimento per up e down
matrice=matrix(c(sum(row.names(up_1_4)%in%id_down_loser),sum(row.names(up_1_4)%in%back_all)-sum(row.names(up_1_4)%in%id_down_loser),sum(id_down_loser%in%back_all)-sum(row.names(up_1_4)%in%id_down_loser),sum(!(back_all%in%row.names(up_1_4))&!(back_all%in%id_down_loser))),nrow=2,ncol=2,byrow = T)
#Fisher test
fisher.test(matrice)
##
## Fisher's Exact Test for Count Data
##
## data: matrice
## p-value = 0.00487
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.6801619 0.9348872
## sample estimates:
## odds ratio
## 0.7979219
#arrichimmento per down e up
matrice=matrix(c(sum(row.names(down_1_4)%in%id_up_loser),sum(row.names(down_1_4)%in%back_all)-sum(row.names(down_1_4)%in%id_up_loser),sum(id_up_loser%in%back_all)-sum(row.names(down_1_4)%in%id_up_loser),sum(!(back_all%in%row.names(down_1_4))&!(back_all%in%id_up_loser))),nrow=2,ncol=2,byrow = T)
#Fisher test
fisher.test(matrice)
##
## Fisher's Exact Test for Count Data
##
## data: matrice
## p-value = 0.006239
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.4867568 0.8964258
## sample estimates:
## odds ratio
## 0.6667054
#di seguito c'e' un esempio copiato da notebook vecchio
#gene ontology analysis
setwd(base_path)
ward_down_no=read.csv("down_ana_last_version.csv")
#go_up_34_function=go_up_34[go_up_34$source=="GO:MF",]
#up_34_function=as.vector(go_up_34_function$term_name)
#go_up_34_process=go_up_34[go_up_34$source=="GO:BP",]
#up_34_process=as.vector(go_up_34_process$term_name)
ward_down_component=ward_down_no[ward_down_no$source=="WP",]
down_component_no=as.vector(ward_down_component$term_name)
down_pvalue_no=as.vector(ward_down_component$adjusted_p_value)
down_output_no=data.frame(down_component_no,down_pvalue_no)
down_cluster_all_function=down_output_no#tabella da vedere
ward_down_component=ward_down_no[ward_down_no$source=="KEGG",]
down_component_no=as.vector(ward_down_component$term_name)
down_pvalue_no=as.vector(ward_down_component$adjusted_p_value)
down_output_no=data.frame(down_component_no,down_pvalue_no)
down_cluster_all_process=down_output_no#tabella da vedere
ward_down_component=ward_down_no[ward_down_no$source=="GO:CC",]
down_component_no=as.vector(ward_down_component$term_name)
down_pvalue_no=as.vector(ward_down_component$adjusted_p_value)
down_output_no=data.frame(down_component_no,down_pvalue_no)
down_vitro_component=down_output_no
#di seguito c'e' un esempio copiato da notebook vecchio
#gene ontology analysis
setwd(base_path)
ward_down_no=read.csv("up_ana_last_version.csv")
#go_up_34_function=go_up_34[go_up_34$source=="GO:MF",]
#up_34_function=as.vector(go_up_34_function$term_name)
#go_up_34_process=go_up_34[go_up_34$source=="GO:BP",]
#up_34_process=as.vector(go_up_34_process$term_name)
ward_down_component=ward_down_no[ward_down_no$source=="TF",]
down_component_no=as.vector(ward_down_component$term_name)
down_pvalue_no=as.vector(ward_down_component$adjusted_p_value)
down_output_no=data.frame(down_component_no,down_pvalue_no)
down_cluster_all_function=down_output_no#tabella da vedere
ward_down_component=ward_down_no[ward_down_no$source=="GO:BP",]
down_component_no=as.vector(ward_down_component$term_name)
down_pvalue_no=as.vector(ward_down_component$adjusted_p_value)
down_output_no=data.frame(down_component_no,down_pvalue_no)
down_cluster_all_process=down_output_no#tabella da vedere
ward_down_component=ward_down_no[ward_down_no$source=="GO:CC",]
down_component_no=as.vector(ward_down_component$term_name)
down_pvalue_no=as.vector(ward_down_component$adjusted_p_value)
down_output_no=data.frame(down_component_no,down_pvalue_no)
up_vitro_component=down_output_no
#di seguito c'e' un esempio copiato da notebook vecchio
#gene ontology analysis
setwd(base_path)
ward_down_no=read.csv("down_vivo_last_version.csv")
#go_up_34_function=go_up_34[go_up_34$source=="GO:MF",]
#up_34_function=as.vector(go_up_34_function$term_name)
#go_up_34_process=go_up_34[go_up_34$source=="GO:BP",]
#up_34_process=as.vector(go_up_34_process$term_name)
ward_down_component=ward_down_no[ward_down_no$source=="GO:MF",]
down_component_no=as.vector(ward_down_component$term_name)
down_pvalue_no=as.vector(ward_down_component$adjusted_p_value)
down_output_no=data.frame(down_component_no,down_pvalue_no)
down_cluster_all_function=down_output_no#tabella da vedere
ward_down_component=ward_down_no[ward_down_no$source=="GO:BP",]
down_component_no=as.vector(ward_down_component$term_name)
down_pvalue_no=as.vector(ward_down_component$adjusted_p_value)
down_output_no=data.frame(down_component_no,down_pvalue_no)
down_cluster_all_process=down_output_no#tabella da vedere
ward_down_component=ward_down_no[ward_down_no$source=="GO:CC",]
down_component_no=as.vector(ward_down_component$term_name)
down_pvalue_no=as.vector(ward_down_component$adjusted_p_value)
down_output_no=data.frame(down_component_no,down_pvalue_no)
down_vivo_component=down_output_no
#di seguito c'e' un esempio copiato da notebook vecchio
#gene ontology analysis
setwd(base_path)
ward_down_no=read.csv("up_vivo_last_version.csv")
#go_up_34_function=go_up_34[go_up_34$source=="GO:MF",]
#up_34_function=as.vector(go_up_34_function$term_name)
#go_up_34_process=go_up_34[go_up_34$source=="GO:BP",]
#up_34_process=as.vector(go_up_34_process$term_name)
ward_down_component=ward_down_no[ward_down_no$source=="GO:MF",]
down_component_no=as.vector(ward_down_component$term_name)
down_pvalue_no=as.vector(ward_down_component$adjusted_p_value)
down_output_no=data.frame(down_component_no,down_pvalue_no)
down_cluster_all_function=down_output_no#tabella da vedere
ward_down_component=ward_down_no[ward_down_no$source=="GO:BP",]
down_component_no=as.vector(ward_down_component$term_name)
down_pvalue_no=as.vector(ward_down_component$adjusted_p_value)
down_output_no=data.frame(down_component_no,down_pvalue_no)
down_cluster_all_process=down_output_no#tabella da vedere
ward_down_component=ward_down_no[ward_down_no$source=="GO:CC",]
down_component_no=as.vector(ward_down_component$term_name)
down_pvalue_no=as.vector(ward_down_component$adjusted_p_value)
down_output_no=data.frame(down_component_no,down_pvalue_no)
up_vivo_component=down_output_no
Extended Figure 10 C
library(UpSetR)
down_vivo_component=down_vivo_component$down_component_no
up_vivo_component=up_vivo_component$down_component_no
down_vitro_component=down_vitro_component$down_component_no
up_vitro_component=up_vitro_component$down_component_no
listInput <- list(up_vitro_component = up_vitro_component, down_vitro_component =down_vitro_component, up_vivo_component = up_vivo_component,down_vivo_component=down_vivo_component)
upset(fromList(listInput), order.by = "freq")

setwd("/Users/gabriele.lubatti/Desktop/Phd/Cell_Competition/Cell_Competition_Github/SC_RNA_Seq_published_dataset/Cheng_et_al_2019")
#bo=read.csv("GSE109071_series_matrix.txt",sep="\t",header = F)
meta_dati_1=read.csv("SraRunTable.txt.csv")
meta_dati=read.table("sample.tsv",sep="\t",header=T)
dati_raw=read.table("GSE109071_read.txt")
age=meta_dati_1$Age
age=(as.vector(age))
age_new=age[age=="embryonic day 6.5"|age=="embryonic day 5.5" |age=="embryonic day 6.25"]
index=age=="embryonic day 6.5"|age=="embryonic day 5.5" |age=="embryonic day 6.25"
sample=meta_dati_1$Sample.Name
sample=as.vector(sample)
accession=meta_dati$Accession
nome=as.vector(meta_dati$Title)
nome=substring(nome,19,nchar(nome))
provo_1=data.frame(accession,nome)
row.names(provo_1)=accession
provo_2=provo_1[sample,]
provo_3=provo_2
dati_raw_1=t(dati_raw)
dati_raw_2=dati_raw_1[as.vector(provo_3$nome)[as.vector(provo_3$nome)%in%colnames(dati_raw)],]
dati_raw_3=t(dati_raw_2)
nome_3=nome_gene$Gene.name[nome_gene$Gene.name%in%row.names(dati_raw_3)]
index_3=nome_gene$Gene.stable.ID[nome_gene$Gene.name%in%row.names(dati_raw_3)]
dati_raw_4=dati_raw_3[nome_3,]
row.names(dati_raw_4)=index_3
nome_index=as.vector(provo_3$nome[index])
dati_raw_5=dati_raw_4[,nome_index]
data_raw_fmk=data_raw[ ,cells_epi_fmk]
data_tpm_fmk=data_tpm[,cells_epi_fmk]
data_raw_dmso=data_raw[,cells_epi_dmso]
data_tpm_dmso=data_tpm[,cells_epi_dmso]
geni_comuni_nuovi=intersect(row.names(dati_raw_5),row.names(data_raw_fmk))
dati_raw_5=dati_raw_5[geni_comuni_nuovi,]
data_raw_fmk=data_raw_fmk[geni_comuni_nuovi,]
data_tpm_fmk=data_tpm_fmk[geni_comuni_nuovi,]
data_raw_dmso=data_raw_dmso[geni_comuni_nuovi,]
data_tpm_dmso=data_tpm_dmso[geni_comuni_nuovi,]
data_pro <- CreateSeuratObject(counts = dati_raw_5, project = "data_wolf")
#data_diego$stim <- "MUC13230"
#data_1_seurat$isnt <- "IPS"
#print(head(x = rownames(x = SIGAH12)))
#print(head(x = colnames(x = SIGAH12)))
#print(SIGAH12)
data_pro<- subset(data_pro, subset = nFeature_RNA > 0)
data_pro <- NormalizeData(data_pro, verbose = FALSE)
data_pro<- FindVariableFeatures(data_pro, selection.method = "vst", nfeatures = 3000)
data_pro<- ScaleData(data_pro, verbose = FALSE)
data_pro <- RunPCA(data_pro, npcs = 30, verbose = FALSE)
data_pro<- RunUMAP(data_pro, reduction = "pca", dims = 1:20,set.seed=42)
## Warning: The following arguments are not used: set.seed
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 12:31:25 UMAP embedding parameters a = 0.9922 b = 1.112
## 12:31:25 Read 1310 rows and found 20 numeric columns
## 12:31:25 Using Annoy for neighbor search, n_neighbors = 30
## 12:31:25 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 12:31:25 Writing NN index file to temp file /var/folders/sh/dxnz33055m57vnkk6rd8d_b561qpqw/T//RtmpkRZ4tm/file110c44240d420
## 12:31:25 Searching Annoy index using 1 thread, search_k = 3000
## 12:31:25 Annoy recall = 100%
## 12:31:26 Commencing smooth kNN distance calibration using 1 thread
## 12:31:27 Initializing from normalized Laplacian + noise
## 12:31:27 Commencing optimization for 500 epochs, with 47972 positive edges
## 12:31:29 Optimization finished
data_pro <- RunTSNE(data_pro, reduction = "pca", dims = 1:20)
data_pro <- FindNeighbors(data_pro, reduction = "pca", dims = 1:20,k.param=20)
## Computing nearest neighbor graph
## Computing SNN
data_pro <- FindClusters(data_pro, resolution =0.003)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 1310
## Number of edges: 34879
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9985
## Number of communities: 3
## Elapsed time: 0 seconds
data_pro$"label"=age_new
library(repr)
library(cowplot)
options(repr.plot.width=8, repr.plot.height=3)
#plot batches and clusters
#p1 <- DimPlot(ips_combined, reduction = "umap", group.by = "stim", dim.1 = 10, dim.2 = 10)
p2 <- DimPlot(data_pro, reduction = "umap", label = TRUE, dims = c(1,2))
p3 <- DimPlot(data_pro, reduction = "pca", label = TRUE, dims = c(1,2))
plt <- plot_grid( p2,p3,rel_widths = c(1.1, 1))
plt

DefaultAssay(data_pro) <- "RNA"
p3 <- DimPlot(data_pro, reduction = "pca", label = TRUE, dims = c(1,2),group.by = "label")
plt <- plot_grid(p3,rel_widths = c(1.1, 1))
plt

p3 <- DimPlot(data_pro, reduction = "umap", label = TRUE, dims = c(1,2),group.by = "label")
p3 <- DimPlot(data_pro, reduction = "umap", label = TRUE, dims = c(1,2),group.by = "label")
plt <- plot_grid(p3,rel_widths = c(1.1, 1))
plt

p3 <- DimPlot(data_pro, reduction = "tsne", label = TRUE, dims = c(1,2))
plt <- plot_grid(p3,rel_widths = c(1.1, 1))
plt

hvg_pro=data_pro@assays$RNA@var.features
DefaultAssay(data_pro) <- "RNA"
PlotGeneExpression <- function(Sobj,gname){
#plot gene expression in UMAP
p1 <- FeaturePlot(Sobj, reduction="pca",features = c(gname))
#plot UMAP colored by cluster
#p2 <- DimPlot(Sobj, reduction = "umap", label = TRUE)
#plt <- plot_grid(p1, p2,rel_widths = c(1.2, 1))
#plot UMAP colored according to IVF/NT condition
p3 <- DimPlot(Sobj, reduction = "pca", group.by = "label")
#p4 <- VlnPlot(Sobj, features = c(gname))
#Arrange them together in a grid
#top_row <- plot_grid(p1, p2, p3)
plot_grid(p1, ncol = 1, nrow=2)
}
# EXAMPLE OF USAGE
#nome_gene$Gene.stable.ID[nome_gene$Gene.name=="Pou5f1"]
gname<-nome_gene$Gene.stable.ID[nome_gene$Gene.name=="Pou5f1"]
PlotGeneExpression(data_pro,gname)

gname<-nome_gene$Gene.stable.ID[nome_gene$Gene.name=="Bmp4"]
PlotGeneExpression(data_pro,gname)

gname<-nome_gene$Gene.stable.ID[nome_gene$Gene.name=="Amn"]
PlotGeneExpression(data_pro,gname)

#Cellule ExE
dati_pca=(Embeddings(data_pro, reduction = "pca")[, 1:2])
dati_pca=as.data.frame(dati_pca)
cellule_exe_1=row.names(dati_pca)[(dati_pca$PC_1 >-5) & (dati_pca$PC_2<-20)]
projec<-as.matrix(GetAssayData(data_pro, slot ="data",assay="RNA"))
logico=projec[nome_gene$Gene.stable.ID[nome_gene$Gene.name=="Bmp4"],]>1.5
cellule_exe_2=colnames(projec)[logico]
cellule_exe=intersect(cellule_exe_1,cellule_exe_2)
#Cellule Epiblasto
dati_pca=(Embeddings(data_pro, reduction = "pca")[, 1:2])
dati_pca=as.data.frame(dati_pca)
dim(dati_pca)
## [1] 1310 2
cellule_epi_2=row.names(dati_pca)[(dati_pca$PC_1>5)&(dati_pca$PC_2<20)]
projec<-as.matrix(GetAssayData(data_pro, slot ="data",assay="RNA"))
logico=projec[nome_gene$Gene.stable.ID[nome_gene$Gene.name=="Pou5f1"],]>1.5
cellule_epi_1=colnames(projec)[logico]
cellule_epi=intersect(cellule_epi_1,cellule_epi_2)
#Cellule VE
dati_pca=(Embeddings(data_pro, reduction = "pca")[, 1:2])
dati_pca=as.data.frame(dati_pca)
dim(dati_pca)
## [1] 1310 2
cellule_ve_2=row.names(dati_pca)[(dati_pca$PC_1<0)&(dati_pca$PC_2<0)]
projec<-as.matrix(GetAssayData(data_pro, slot ="data",assay="RNA"))
logico=projec[nome_gene$Gene.stable.ID[nome_gene$Gene.name=="Amn"],]>1.5
cellule_ve_1=colnames(projec)[logico]
cellule_ve=intersect(cellule_ve_1,cellule_ve_2)
cellule_all=c(cellule_epi,cellule_exe,cellule_ve)
cellule_all=cellule_all[isUnique(cellule_all)]
logico_2=as.vector(provo_3$nome)%in%cellule_all
age_spero=age[logico_2]
cell_all_new=as.vector(provo_3$nome)[as.vector(provo_3$nome)%in%cellule_all]
#lavoro con cellule_epi_new ,age_spero e projec_new
projec_new=projec[,cell_all_new]
dati_raw_6=dati_raw_5[,cell_all_new]
cell_type=rep(0,length(cell_all_new))
cell_type[cell_all_new%in%cellule_epi]="Epiblast"
cell_type[cell_all_new%in%cellule_exe]="ExE"
cell_type[cell_all_new%in%cellule_ve]="VE"
Epiblast cells at stages E5.5, E6.25 and E6.5 (Cheng et al)
dati_pca=(Embeddings(data_pro, reduction = "pca")[, 1:2])
dati_pca=as.data.frame(dati_pca)
dim(dati_pca)
## [1] 1310 2
cellule_epi_2=row.names(dati_pca)[dati_pca$PC_1>5&dati_pca$PC_2<20]
projec<-as.matrix(GetAssayData(data_pro, slot ="data",assay="RNA"))
logico=projec[nome_gene$Gene.stable.ID[nome_gene$Gene.name=="Pou5f1"],]>1.5
cellule_epi_1=colnames(projec)[logico]
cellule_epi=intersect(cellule_epi_1,cellule_epi_2)
logico_2=as.vector(provo_3$nome)%in%cellule_epi
age_spero=age[logico_2]
cell_epi_new=as.vector(provo_3$nome)[as.vector(provo_3$nome)%in%cellule_epi]
#lavoro con cellule_epi_new ,age_spero e projec_new
projec_new=projec[,cell_epi_new]
dati_raw_6=dati_raw_5[,cell_epi_new]
data_de=dati_raw_6
#tengo solo i geni con media tpm maggiore di 10 nei miei vecchi dati
geni=row.names(data_de)[row.names(data_de)%in%row.names(data_tpm)]
dati_confronto=data_de[geni,]
setwd(base_path)
#lunghezze=getGeneLengthAndGCContent(row.names(dati_confronto), org="mmu")
#lunghezze_data=as.data.frame(lunghezze)
#save(lunghezze_data,file="lunghezze_data.Rda")
load("lunghezze_data.Rda")
#lunghezze_data_wolf=lunghezze_data[geni,]
geni_nuovo_nuovo=row.names(lunghezze_data)[!(is.na(lunghezze_data$length))]
lunghezze_vere=lunghezze_data$length[!(is.na(lunghezze_data$length))]
geni_tengo=geni_nuovo_nuovo[geni_nuovo_nuovo%in%row.names(dati_confronto)]
lunghezze_vere_vere=lunghezze_vere[geni_nuovo_nuovo%in%row.names(dati_confronto)]
dati_confronto_nuovi=dati_confronto[geni_tengo,]
sce_jon<-SingleCellExperiment(list(counts=as.matrix(dati_confronto_nuovi)))
tpm_proiezione_jon=calculateTPM(sce_jon,lengths=lunghezze_vere_vere)
sce_pro_jon<- SingleCellExperiment(assays = list(normcounts = as.matrix(tpm_proiezione_jon)))
logcounts(sce_pro_jon) <- log10(normcounts(sce_pro_jon) + 1)
rowData(sce_pro_jon)$feature_symbol <- rownames(sce_pro_jon)
data_6_5 <- CreateSeuratObject(counts = data_de, project = "data_wolf")
#data_diego$stim <- "MUC13230"
#data_1_seurat$isnt <- "IPS"
#print(head(x = rownames(x = SIGAH12)))
#print(head(x = colnames(x = SIGAH12)))
#print(SIGAH12)
data_6_5<- subset(data_6_5, subset = nFeature_RNA > 0)
data_6_5 <- NormalizeData(data_6_5, verbose = FALSE)
data_6_5<- FindVariableFeatures(data_6_5, selection.method = "vst", nfeatures = 1000)
HVGC=VariableFeatures(data_6_5)
dist_computation=function(a,b){
log_data_tpm <- log10(a+1)[b,]
#save(log_data_tpm,file="chosen_exprs1_novembre.Rda")
#load(file="chosen_exprs1_novembre.Rda")
log_data_tpm<-as.matrix(log_data_tpm)
dissimilarity <- sqrt((1-cor((log_data_tpm), method = "spearman"))/2)
#save(dissimilarity,file="dissimilarity_all.Rda")
#Umap by medium
return(dissimilarity)
}
geni_comuni=HVGC[HVGC%in%row.names(tpm_proiezione_jon)]
tpm_proiezione_jon=tpm_proiezione_jon[geni_comuni,]
data_tpm_fmk_jon=data_tpm_fmk[geni_comuni,]
all_data=cbind(data_tpm_fmk_jon,tpm_proiezione_jon)
#all_data=data_tpm_fmk
#load("HVGC_giusto_no_un.Rda")
#load("dissimilarity_all.Rda")
#dissimilarity_all=dissimilarity
#log_data_tpm <- log10(data_tpm[HVGC,]+1)
#log_data_tpm<-as.matrix(log_data_tpm)
#dissimilarity_all <- sqrt((1-cor((log_data_tpm), method = "spearman"))/2)
dissimilarity_all=dist_computation(all_data,geni_comuni)
#dissimilarity_fmk=dist_computation(fmk,hvg_wolf)
Extended Figure 3 E and F
set.seed(1000)
dissimilarity_fmk=dissimilarity_all[colnames(data_tpm_fmk_jon),colnames(tpm_proiezione_jon)]
dissimilarity_jon=dissimilarity_all[colnames(tpm_proiezione_jon),colnames(tpm_proiezione_jon)]
dm=DiffusionMap(dissimilarity_jon,distance='cosine')
## Warning in DiffusionMap(dissimilarity_jon, distance = "cosine"): You have 520
## genes. Consider passing e.g. n_pcs = 50 to speed up computation.
#fmk=fmk[hvg_fmk,]
#dm=DiffusionMap(t(fmk))
DC1=dm$DC1
DC2=dm$DC2
#Cluster=as.factor(Cluster_col)
#clusters_fmk=as.vector(wolf_6_5$label)
#clusters=meta_dati_course[cells_epi_fmk,]$cluster
clusters=age_spero
clusters[age_spero=="embryonic day 5.5"]="5"
clusters[age_spero=="embryonic day 6.25" ]="6"
clusters[age_spero=="embryonic day 6.5" ]="7"
#clusters[clusters==1]="Normal (winner) Epiblast"
#clusters[clusters==3]="Intermediate"
#clusters[clusters==4]="Loser Epiblast"
Cluster=as.factor(clusters)
tre=data.frame(DC1,DC2,Cluster)
#dpt <- DPT(dm,tips=20)
#dpt <- DPT(dm)
#pdf("figure_paper_pdf/diffusion_map_all.pdf",width=8,height=5)
diffu<-ggplot(tre, aes(x=DC1 , y=DC2 , color= Cluster)) + geom_point()+
labs(x = "DC1",y="DC2",col="Clusters") + ggtitle("Diffusion Map")+
scale_color_manual(breaks = c("5","6","7")
,labels = c("Epiblast 5.5"
,"Epiblast 6.25"
,"Epiblast 6.5"
)
,values=c("grey","yellow","red" ))+
labs(col="Cluster")
diffu

#dpt <- DPT(dm)
dpt <- DPT(dm,tips=150)
#save(dpt,file="dpt.Rda")
#str(dpt)
tips(dpt)
## EB_2223 EB_950 EB_2054
## 150 333 102
plot(dpt, col_by = 'DPT333')

#wolf<-as.matrix(GetAssayData(wolf_6_5, slot ="data",assay="RNA"))
#wolf=wolf[hvg_fmk,]
pred=dm_predict(dm,(dissimilarity_fmk))
prima=c(dm$DC1,pred[,1])
seconda=c(dm$DC2,pred[,2])
clusters_pro_fmk=meta_dati_course[cells_epi_fmk,]$cluster
#clusters_pro=meta_dati_course[cells_epi_dmso,]$cluster
#clusters_pro[clusters_pro==1]="5"
#clusters_pro[clusters_pro==3]="6"
#clusters_pro=rep(5,length(colnames(data_tpm_dmso)))
#Cluster_pro=as.factor(clusters_pro)
colore_plot=c(clusters,clusters_pro_fmk)
colore_plot=as.factor(colore_plot)
dati_plot=data.frame(prima,seconda,colore_plot)
diffu<-ggplot(dati_plot, aes(x=prima , y=seconda , color=(colore_plot))) + geom_point()+labs(x = "DC1",y="DC2") + ggtitle("Diffusion Map")+labs(col="Cluster")+scale_color_manual(breaks = c("1","3","4","5","6","7")
,labels = c("Normal (winner) Epiblast"
,"Intermediate"
,"Loser Epiblast "
,"Epiblast 5.5"
,"Epiblast 6.25"
,"Epiblast 6.5"
)
,values=c("#DD6400","#0000FF","#006400","grey","yellow","red" ))
diffu

time_dati_nuovi=dpt$DPT333
projection_dist_ten=function (dm, new_dcs = NULL, ..., new_data, verbose = FALSE)
{
if (is.null(new_dcs))
new_dcs <- dm_predict(dm, new_data, ..., verbose = verbose)
else if (!missing(new_data))
stop("only pass one of new_dcs and new_data")
evs <- eigenvectors(dm)
nns <- find_knn(evs, 10, query = as.matrix(new_dcs))
nn_dist <- nns$index[, 1:10]
#nn_dist
}
dist=projection_dist_ten(dm,new_dcs=pred)
k10=matrix(0,ncol=10,nrow = length(colnames(data_tpm_fmk_jon)))
for( i in 1:length(colnames(data_tpm_fmk_jon))){
k10[i,]=time_dati_nuovi[dist[i,]]
#print(i)
}
time_projec_fmk=apply(k10,1,mean)
dfForPlot <- data.frame(clusters_fmk = as.character(c(clusters,clusters_pro_fmk))
,exp=c(time_dati_nuovi,time_projec_fmk)
,Cluster = as.factor(as.character(c(clusters,clusters_pro_fmk))))
ggplot(data =dfForPlot
,aes(x = clusters_fmk
,y = exp
,col = Cluster )) +
geom_boxplot(width=0.9,alpha=0.8) +
scale_color_manual(breaks = c("1","3","4","5","6","7")
,labels = c("Normal (winner) Epiblast"
,"Intermediate"
,"Loser Epiblast "
,"Epiblast 5.5"
,"Epiblast 6.25"
,"Epiblast 6.5"
)
,values=c("#DD6400","#0000FF","#006400" ,"grey","yellow","red")) +
ylab("DPT")+xlab(NULL) +
geom_sina()+
ggtitle("Pseudo time FMK cells (274) and epiblast cells E5.5, E6.25 and E6.5 (520)")+ylim(c(0,max(time_dati_nuovi)))+theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank())

#pdf("losing_score_dmso_projection.pdf",useDingbats = F)
Projection DMSO epiblast cells
geni_comuni=HVGC[HVGC%in%row.names(tpm_proiezione_jon)]
tpm_proiezione_jon=tpm_proiezione_jon[geni_comuni,]
data_tpm_dmso_jon=data_tpm_dmso[geni_comuni,]
all_data=cbind(data_tpm_dmso_jon,tpm_proiezione_jon)
#all_data=data_tpm_fmk
#load("HVGC_giusto_no_un.Rda")
#load("dissimilarity_all.Rda")
#dissimilarity_all=dissimilarity
#log_data_tpm <- log10(data_tpm[HVGC,]+1)
#log_data_tpm<-as.matrix(log_data_tpm)
#dissimilarity_all <- sqrt((1-cor((log_data_tpm), method = "spearman"))/2)
dissimilarity_all=dist_computation(all_data,geni_comuni)
#dissimilarity_fmk=dist_computation(fmk,hvg_wolf)
set.seed(1001)
dissimilarity_dmso=dissimilarity_all[colnames(data_tpm_dmso_jon),colnames(tpm_proiezione_jon)]
dissimilarity_jon=dissimilarity_all[colnames(tpm_proiezione_jon),colnames(tpm_proiezione_jon)]
dm=DiffusionMap(dissimilarity_jon,distance="cosine")
## Warning in DiffusionMap(dissimilarity_jon, distance = "cosine"): You have 520
## genes. Consider passing e.g. n_pcs = 50 to speed up computation.
#dm=DiffusionMap(dissimilarity_jon)
#fmk=fmk[hvg_fmk,]
#dm=DiffusionMap(t(fmk))
DC1=dm$DC1
DC2=dm$DC2
#Cluster=as.factor(Cluster_col)
#clusters_fmk=as.vector(wolf_6_5$label)
#clusters=meta_dati_course[cells_epi_fmk,]$cluster
clusters=age_spero
clusters[age_spero=="embryonic day 5.5"]="5"
clusters[age_spero=="embryonic day 6.25" ]="6"
clusters[age_spero=="embryonic day 6.5" ]="7"
#clusters[clusters==1]="Normal (winner) Epiblast"
#clusters[clusters==3]="Intermediate"
#clusters[clusters==4]="Loser Epiblast"
Cluster=as.factor(clusters)
tre=data.frame(DC1,DC2,Cluster)
#dpt <- DPT(dm,tips=20)
#dpt <- DPT(dm)
#pdf("figure_paper_pdf/diffusion_map_all.pdf",width=8,height=5)
diffu<-ggplot(tre, aes(x=DC1 , y=DC2 , color= Cluster)) + geom_point()+
labs(x = "DC1",y="DC2",col="Clusters") + ggtitle("Diffusion Map")+
scale_color_manual(breaks = c("5","6","7")
,labels = c("Epiblast 5.5"
,"Epiblast 6.25"
,"Epiblast 6.5"
)
,values=c("grey","yellow","red" ))+
labs(col="Cluster")
diffu

#dpt <- DPT(dm)
dpt <- DPT(dm,tips=150)
#save(dpt,file="dpt.Rda")
#str(dpt)
tips(dpt)
## EB_2223 EB_950 EB_2054
## 150 333 102
plot(dpt, col_by = 'DPT333')

#wolf<-as.matrix(GetAssayData(wolf_6_5, slot ="data",assay="RNA"))
#wolf=wolf[hvg_fmk,]
pred=dm_predict(dm,(dissimilarity_dmso))
prima=c(dm$DC1,pred[,1])
seconda=c(dm$DC2,pred[,2])
clusters_pro_dmso=meta_dati_course[cells_epi_dmso,]$cluster
#clusters_pro=meta_dati_course[cells_epi_dmso,]$cluster
#clusters_pro[clusters_pro==1]="5"
#clusters_pro[clusters_pro==3]="6"
#clusters_pro=rep(5,length(colnames(data_tpm_dmso)))
#Cluster_pro=as.factor(clusters_pro)
colore_plot=c(clusters,clusters_pro_dmso)
colore_plot=as.factor(colore_plot)
dati_plot=data.frame(prima,seconda,colore_plot)
diffu<-ggplot(dati_plot, aes(x=prima , y=seconda , color=(colore_plot))) + geom_point()+labs(x = "DC1",y="DC2") + ggtitle("Diffusion Map")+labs(col="Cluster")+scale_color_manual(breaks = c("1","3","5","6","7")
,labels = c("Normal (winner) Epiblast"
,"Intermediate"
,"Epiblast 5.5"
,"Epiblast 6.25"
,"Epiblast 6.5"
)
,values=c("#DD6400","#0000FF","grey","yellow","red" ))
diffu

time_dati_nuovi=dpt$DPT333
projection_dist_ten=function (dm, new_dcs = NULL, ..., new_data, verbose = FALSE)
{
if (is.null(new_dcs))
new_dcs <- dm_predict(dm, new_data, ..., verbose = verbose)
else if (!missing(new_data))
stop("only pass one of new_dcs and new_data")
evs <- eigenvectors(dm)
nns <- find_knn(evs, 10, query = as.matrix(new_dcs))
nn_dist <- nns$index[, 1:10]
#nn_dist
}
dist=projection_dist_ten(dm,new_dcs=pred)
k10=matrix(0,ncol=10,nrow = length(colnames(data_tpm_dmso_jon)))
for( i in 1:length(colnames(data_tpm_dmso_jon))){
k10[i,]=time_dati_nuovi[dist[i,]]
#print(i)
}
time_projec_dmso=apply(k10,1,mean)
Extended Figure 3 G
dfForPlot <- data.frame(clusters_tutto = as.character(c(clusters,clusters_pro_fmk,clusters_pro_dmso))
,exp=c(time_dati_nuovi,time_projec_fmk,time_projec_dmso)
,Cluster = as.factor(as.character(c(clusters,clusters_pro_fmk,clusters_pro_dmso))))
clusters_dmso=clusters_pro_dmso
clusters_fmk=clusters_pro_fmk
clusters_ori=clusters
clusters_dmso[clusters_pro_dmso%in%3]="5"
clusters_dmso[clusters_pro_dmso%in%1]="4"
clusters_fmk[clusters_pro_fmk%in%4]="3"
clusters_fmk[clusters_pro_fmk%in%3]="2"
clusters_fmk[clusters_pro_fmk%in%1]="1"
clusters_ori[clusters%in%5]="6"
clusters_ori[clusters%in%6]="7"
clusters_ori[clusters%in%7]="8"
clusters_all=c(clusters_ori,clusters_fmk,clusters_dmso)
nuovo_ordine <- factor(clusters_all,
levels = c("1","2","3","4","5","6","7","8",ordered = TRUE))
#all_clu=c(ci_name,dmso_name)
ggplot(data =dfForPlot
,aes(x = nuovo_ordine
,y = dfForPlot$exp
,col = nuovo_ordine )) +
geom_boxplot(width=0.9,alpha=0.8) +
scale_color_manual(breaks = c("1","2","3","4","5","6","7","8")
,labels = c("Winner Epiblast-CI"
,"Intermediate-CI"
,"Loser Epiblast-CI"
,"Winner Epiblast-DMSO"
,"Intermediate-DMSO"
,"Epiblast 5.5"
,"Epiblast 6.25"
,"Epiblast 6.5"
)
,values=c("#DD6400","#0000FF","#006400","#DD6400","#0000FF", "grey","yellow","red")) +
ylab("DPT")+xlab(NULL) + labs(col="Cluster")+
geom_sina()+
ggtitle("Pseudo time CI (274), DMSO (238) and epiblast cells E5.5, E6.25 and E6.5 (520)")+ylim(c(0,max(time_dati_nuovi)))+
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank())
## Warning: Use of `dfForPlot$exp` is discouraged. Use `exp` instead.
## Warning: Use of `dfForPlot$exp` is discouraged. Use `exp` instead.

#pdf("losing_score_dmso_projection.pdf",useDingbats = F)
Extended Figure 2 A
data_raw_fmk=data_raw[,cells_epi_fmk]
data_raw_dmso=data_raw[,cells_epi_dmso]
meta_dati_course_fmk=meta_dati_course[cells_epi_fmk,]
meta_dati_course_dmso=meta_dati_course[cells_epi_dmso,]
mm.pairs <- readRDS(system.file("exdata", "mouse_cycle_markers.rds", package="scran"))
out1<-list(g1=mm.pairs$G1,g2=mm.pairs$G2M,s=mm.pairs$S)
cyc<-cyclone(as.matrix(data_raw_fmk),out1)
score_r<-cyc$scores
g1_max<-score_r$g1>= 0.5 & score_r$g1>=score_r$g2
g2_max<-score_r$g2>=0.5 & score_r$g2>score_r$g1
s_max<-score_r$g1<0.5 & score_r$g2<0.5
sequenza_fmk<-rep(0,length(cells_epi_fmk))
sequenza_fmk[g1_max]="G1"
sequenza_fmk[s_max]="S"
sequenza_fmk[g2_max]="G2M"
provo_fmk=as.factor(sequenza_fmk)
provo_fmk <- factor(provo_fmk, levels = c("G2M", "S", "G1"))
cyc<-cyclone(as.matrix(data_raw_dmso),out1)
score_r<-cyc$scores
g1_max<-score_r$g1>= 0.5 & score_r$g1>=score_r$g2
g2_max<-score_r$g2>=0.5 & score_r$g2>score_r$g1
s_max<-score_r$g1<0.5 & score_r$g2<0.5
sequenza_dmso<-rep(0,length(cells_epi_dmso))
sequenza_dmso[g1_max]="G1"
sequenza_dmso[s_max]="S"
sequenza_dmso[g2_max]="G2M"
provo_dmso=as.factor(sequenza_dmso)
provo_dmso <- factor(provo_dmso, levels = c("G2M", "S", "G1"))
provo_all=c(sequenza_fmk,sequenza_dmso)
provo_all <- factor(provo_all, levels = c("G2M", "S", "G1"))
ci_name=meta_dati_course_fmk$cluster
ci_name[ci_name==3]="Intermediate-CI"
ci_name[ci_name==4]="Los Epi-CI"
ci_name[ci_name==1]="Win Epi-CI"
dmso_name=meta_dati_course_dmso$cluster
dmso_name[dmso_name==3]="Intermediate-DMSO"
dmso_name[dmso_name==4]="Los Epi-DMSO"
dmso_name[dmso_name==1]="Win Epi-DMSO"
name_full=c(ci_name,dmso_name)
nuovo_ordine <- factor(name_full,
levels = c("Win Epi-CI","Win Epi-DMSO","Intermediate-CI","Intermediate-DMSO","Los Epi-CI","Los Epi-DMSO",ordered = TRUE))
#all_clu=c(ci_name,dmso_name)
dfForPlot <- data.frame(Phases =provo_all
,clu =nuovo_ordine)
pnm <-ggplot(dfForPlot, aes(y=Phases,x=clu))
pnm +geom_bar( aes(fill = Phases,y = (..count..)/sum(..count..)),position="fill")+ ggtitle(paste( as.character("Cell cycle CI and DMSO cells")))+ylab("Fraction of cells")+xlab("Cluster")

dfForPlot <- data.frame(Phases =sequenza_fmk
,clu =as.factor(meta_dati_course_fmk$cluster))
dfForPlot <- data.frame(Phases =provo_all
,clu =nuovo_ordine)
pnm <-ggplot(dfForPlot, aes(y=Phases,x=clu))
pnm +geom_bar( aes(fill = Phases,y = (..count..)))+ ggtitle(paste( as.character("Cell cycle CI and DMSO cells")))+ylab("Number of cells")+xlab("Cluster")

Heteroplasmy
setwd(base_path)
load("meta_dati.Rda")
nomi=read.csv("old_2_cell_names.csv")
nomi=nomi$X0
nomi=as.vector(nomi)
setwd(base_path)
all_pos=read.csv("old_2_mt_all.csv")
all_pos=all_pos[,-1]
colnames(all_pos)=substring(colnames(all_pos),2,nchar(colnames(all_pos)))
non_so=gsub("[.]", "_", colnames(all_pos))
colnames(all_pos)=non_so
row.names(all_pos)=substring(nomi,12,nchar(nomi)-25)
vedere=substring(colnames(all_pos),1,nchar(colnames(all_pos))-3)
vedere_numb=substring(vedere,1,nchar(vedere)-1)
#nomi cellule e cluster
setwd(base_path)
load("provo_ward_no_un.Rda")
load("nomi_full.Rda")#name of the old bam file
#load("nomi_nuovi.Rda")#name on array express
#load("nomi_vecchi.Rda")#
nomi_nuovi=meta_dati$nomi_nuovi
nomi_vecchi=meta_dati$nomi_vecchi
row.names(provo_ward_no_un)=provo_ward_no_un$Cellname
provo_ward_no_un=provo_ward_no_un[nomi_full,]
load("provo.Rda")
cells_epi=row.names(provo)
row.names(provo_ward_no_un)=cells_epi
provo_1=provo_ward_no_un[!(grepl("2",provo_ward_no_un$cluster)),]
provo_1$nomi_vecchi=row.names(provo_1)
old=nomi_vecchi[nomi_vecchi%in%row.names(provo_1)]
provo_2=provo_1[old,]
new=nomi_nuovi[nomi_vecchi%in%row.names(provo_1)]
provo_2$nomi_nuovi=new
row.names(provo_2)=provo_2$Cellname
cell_common=row.names(all_pos)[row.names(all_pos)%in%row.names(provo_2)]
provo_3=provo_2[cell_common,]
provo_4=provo_3[grep("FMK",provo_3$condizione),]
cells_fmk=row.names(provo_4)
provo_5=provo_3[grep("DMSO",provo_3$condizione),]
cells_dmso=row.names(provo_5)
cell_clu_1=row.names(provo_4)[provo_4$cluster==1]
cell_clu_3=row.names(provo_4)[provo_4$cluster==3]
cell_clu_4=row.names(provo_4)[provo_4$cluster==4]
all_posizioni=all_pos
vedere_all=non_so
vedere=vedere_numb
entropy_method=function(all_pos,cell_stage,vedere_all,vedere_numb,number_reads,number_positions,filtering=1,selection=FALSE,cell_clu_1=NULL,cell_clu_3=NULL,cell_clu_4=NULL){
all_posizioni=all_pos[cell_stage,]
vedere_all=non_so
vedere=vedere_numb
all_pos=all_pos[cell_stage,]
prova=matrix(0,ncol=length(unique(vedere)),nrow = length(row.names(all_posizioni)))
colnames(prova)=unique(vedere)
row.names(prova)=row.names(all_posizioni)
colnames(all_posizioni)=vedere
for (i in unique(vedere)){
prova[,i]=apply(all_posizioni[,colnames(all_posizioni)%in%i],1,function(x){
x=as.numeric(as.vector(x))
return(sum(x))
})
}
somma=prova
#save(somma,file="somma_antonio.Rda")
#load(file="somma_antonio.Rda")
#First filtering step
somma=as.matrix(somma)
pos_cov=apply(somma,1,function(x){sum(x>number_reads)})
dfForPlot <- data.frame(clusters =row.names(somma)
,genes_expression = as.numeric(pos_cov)
)
p <- ggplot(data =dfForPlot
,aes(x = clusters
,y = genes_expression
))
p+geom_point()+
labs(col="Cluster")+labs(x="Cell")+labs(y="Number of position covered") +ggtitle("QC")
pos_cov_20000=pos_cov[pos_cov>number_positions]
somma_qc=somma[pos_cov>number_positions,]
##Second filtering step
if (filtering==1){
all_pos_update=apply(somma_qc,2,function(x){
if(sum(x>number_reads)==length(row.names(somma_qc))){return (TRUE)}
else{return (FALSE)}
})
somma_qc_qc=somma_qc[,all_pos_update]}
if (filtering==2){
all_pos_update=apply(somma_qc,2,function(x){
clu_1=x[row.names(somma_qc)%in%cell_clu_1]
clu_3=x[row.names(somma_qc)%in%cell_clu_3]
clu_4=x[row.names(somma_qc)%in%cell_clu_4]
if((sum(clu_1>number_reads)>0.5*length(clu_1))&(sum(clu_3>number_reads)>0.5*length(clu_3))&(sum(clu_4>number_reads)>0.5*length(clu_4))){return (TRUE)}
else{return (FALSE)}
})
somma_qc_qc=somma_qc[,all_pos_update]
}
if(selection==TRUE){
test_cell=apply(somma_qc_qc,1,function(x){sum(x>number_reads)})
cell_common=row.names(somma_qc_qc)[which(test_cell==dim(somma_qc_qc)[2])]
somma_qc_qc=somma_qc_qc[cell_common,]
}
indici=apply(somma_qc_qc,2,function(x){
logico=which(x>number_reads)
return((logico))
})
#save(somma_qc_qc,file="somma_qc_qc_antonio.Rda")
#normalization
all_posizioni_qc=all_posizioni[row.names(somma_qc_qc),colnames(all_posizioni)%in%colnames(somma_qc_qc)]
#save(all_posizioni_qc,file="all_posizioni_qc_antonio.Rda")
prova_new=matrix(0,ncol=length(colnames(all_posizioni_qc)),nrow = length(row.names(all_posizioni_qc)))
colnames(prova_new)=vedere_all[colnames(all_posizioni)%in%colnames(somma_qc_qc)]
row.names(prova_new)=row.names(somma_qc_qc)
colnames(all_posizioni_qc)=vedere_all[colnames(all_posizioni)%in%colnames(somma_qc_qc)]
for (i in colnames(all_posizioni_qc)){
prova_new[,i]=all_posizioni_qc[,i]/(somma_qc_qc[,substring(i,1,nchar(i)-4)]+0.0000000001)
}
all_posizioni_qc_norm=prova_new
###Compute Entropy
prova_5=matrix(0,ncol=length(colnames(somma_qc_qc)),nrow=length(row.names(all_posizioni_qc_norm)))
colnames(prova_5)=colnames(somma_qc_qc)
row.names(prova_5)=row.names(all_posizioni_qc_norm)
colnames(all_posizioni_qc_norm)=vedere[vedere_all%in%colnames(all_posizioni_qc)]
for (i in colnames(somma_qc_qc)){
prova_5[,i]=apply(all_posizioni_qc_norm[,colnames(all_posizioni_qc_norm)%in%i],1,function(x){
x=as.numeric(as.vector(x))
if (sum(x)==0){return(-1)}
else{
return(1-max(x))
}
})
}
prova_entro=prova_5
entropia=prova_entro
return(list(entropia,prova_new,indici))
}
plot_heteroplasmy_cell_competition=function(position,entropia_matrix,cluster,indici){
position=colnames(entropia_matrix)[colnames(entropia_matrix)==position]
cluster=cluster[as.numeric(indici[[position]])]
mutazione=entropia_matrix[as.numeric(indici[[position]]),position]
dati_scatter_1=data.frame(as.character(cluster),mutazione)
ggplot(data =dati_scatter_1
,aes(x =factor(cluster,levels=c("Winner Epiblast","Intermediate","Loser Epiblast"),ordered = TRUE)
,y = mutazione
,col =factor(cluster,levels=c("Winner Epiblast","Intermediate","Loser Epiblast"),ordered = TRUE))) +
geom_boxplot(width=0.9,alpha=0.8)+
labs(col="Cluster")+labs(x="Cluster")+labs(y="Heteroplasmy") +ggtitle(position) +scale_color_manual(values=c("#DD6400", "#0000FF","#006400"),labels=c("Winner Epiblast","Intermediate", "Loser Epiblast"))
}
plot_allele_frequency_cell_competition=function(position,entropia_matrix,allele_matrix,cluster,indici,number){
position=colnames(entropia_matrix)[colnames(entropia_matrix)==position]
cluster=cluster[as.numeric(indici[[position]])]
a=position
allele=colnames(allele_matrix)[grep(a,colnames(allele_matrix))]
mutazione_1=allele_matrix[as.numeric(indici[[position]]),allele[1]]
dati_scatter_1=data.frame(cluster,mutazione_1)
plot_1=ggplot(data =dati_scatter_1
,aes(x =factor(cluster,levels=c("Winner Epiblast","Intermediate","Loser Epiblast"),ordered = TRUE)
,y = mutazione_1
,col =factor(cluster,levels=c("Winner Epiblast","Intermediate","Loser Epiblast"),ordered = TRUE))) +
geom_boxplot(width=0.9,alpha=0.8) +
geom_point()+labs(col="Cluster")+labs(y="Allele Frequency")+labs(col="Cluster")+labs(x="Cluster")+labs(y="Allele frequency") +ggtitle(allele[1])+ylim(0,1)+scale_color_manual(values=c("#DD6400", "#0000FF","#006400"),labels=c("Winner Epiblast","Intermediate", "Loser Epiblast"))+theme(axis.text=element_text(size=number))
mutazione_2=allele_matrix[as.numeric(indici[[position]]),allele[2]]
dati_scatter_2=data.frame(cluster,mutazione_2)
plot_2=ggplot(data =dati_scatter_2
,aes(x =factor(cluster,levels=c("Winner Epiblast","Intermediate","Loser Epiblast"),ordered = TRUE)
,y = mutazione_2
,col =factor(cluster,levels=c("Winner Epiblast","Intermediate","Loser Epiblast"),ordered = TRUE))) +
geom_boxplot(width=0.9,alpha=0.8) +
geom_point()+labs(col="Cluster")+labs(y="Allele Frequency")+labs(col="Cluster")+labs(x="Cluster")+labs(y="Allele frequency") +ggtitle(allele[2])+ylim(0,1)+scale_color_manual(values=c("#DD6400", "#0000FF","#006400"),labels=c("Winner Epiblast","Intermediate", "Loser Epiblast"))+theme(axis.text=element_text(size=number))
mutazione_3=allele_matrix[as.numeric(indici[[position]]),allele[3]]
dati_scatter_3=data.frame(cluster,mutazione_3)
plot_3=ggplot(data =dati_scatter_3
,aes(x =factor(cluster,levels=c("Winner Epiblast","Intermediate","Loser Epiblast"),ordered = TRUE)
,y = mutazione_3
,col =factor(cluster,levels=c("Winner Epiblast","Intermediate","Loser Epiblast"),ordered = TRUE))) +
geom_boxplot(width=0.9,alpha=0.8) +
geom_point()+labs(col="Cluster")+labs(y="Allele Frequency")+labs(col="Cluster")+labs(x="Cluster")+labs(y="Allele frequency") +ggtitle(allele[3])+ylim(0,1)+scale_color_manual(values=c("#DD6400", "#0000FF","#006400"),labels=c("Winner Epiblast","Intermediate", "Loser Epiblast"))+theme(axis.text=element_text(size=number))
mutazione_4=allele_matrix[as.numeric(indici[[position]]),allele[4]]
dati_scatter_4=data.frame(cluster,mutazione_4)
plot_4=ggplot(data =dati_scatter_4
,aes(x =factor(cluster,levels=c("Winner Epiblast","Intermediate","Loser Epiblast"),ordered = TRUE)
,y = mutazione_4
,col =factor(cluster,levels=c("Winner Epiblast","Intermediate","Loser Epiblast"),ordered = TRUE))) +
geom_boxplot(width=0.9,alpha=0.8) +
geom_point()+labs(col="Cluster")+labs(y="Allele Frequency")+labs(col="Cluster")+labs(x="Embryo")+labs(y="Allele frequency") +ggtitle(allele[4])+ylim(0,1)+scale_color_manual(values=c("#DD6400", "#0000FF","#006400"),labels=c("Winner Epiblast","Intermediate", "Loser Epiblast"))+theme(axis.text=element_text(size=number))
grid.arrange(plot_1,plot_2,plot_3,plot_4 ,
ncol = 2, nrow = 2)
}
gam_fit_cell_competition=function(entropia_matrix,indici,time,min_heteroplasmy=0.01,min_cells=10){
position=colnames(entropia_matrix)
considero=rep(0,length(position))
for (i in 1:length(position)){
Y=data.frame(t(entropia_matrix[as.numeric(indici[[position[i]]]),position[i]]))
considero[i]=sum(Y>min_heteroplasmy)
considero[considero]
}
entropia_matrix_select=entropia_matrix[,which(considero>min_cells)]
position=colnames(entropia_matrix_select)
p_value=rep(0,length(position))
for (i in 1:length(position)){
Y=data.frame(t(entropia_matrix_select[as.numeric(indici[[position[i]]]),position[i]]))
t=time[as.numeric(indici[[position[i]]])]
#Run GAM with LOESS function for all the genes: save the p-value of parametric ANOVA
#and the fitted expression values for each genes
gam.res <- apply(Y, 1, function(z){
z=as.numeric(as.vector(z))
d <- data.frame(z=z, t=t)
tmp <- gam(z ~ lo(t), data=d)
p <- summary(tmp)[4][[1]][1,5] #p-value parametric ANOVA
f<-fitted(tmp)
#p2 <- summary(tmp)[3][[1]][2,3] #p-value non-parametric ANOVA
c(p,f)
#c(p1,p2) #If we want to return both the parametric and non-parametric p-values
})
genes.table<-data.frame(genes.names=rownames(Y))
genes.table$pvals<-gam.res[1,]
genes.table$FDR<-p.adjust(genes.table$pvals, method="fdr") #Adjust p-values
genes.table<-genes.table[order(genes.table$FDR),] #Order the genes according to the FDR
genes.table$genes.names <- as.character(genes.table$genes.names)
row.names(genes.table)=genes.table$genes.names
results.gam.tot=genes.table
#Save the matrix of fitted values
gam.fitted<-gam.res[-1,]
#results_gam=results.gam.tot['genes.names']
results_gam_nomi=results.gam.tot$genes.names
gam_fitted=gam.fitted
p_value[i]=results.gam.tot$pvals
}
fdr_update=p.adjust(p_value, method="fdr")
all_fdr=data.frame(fdr_update,colnames(entropia_matrix_select))
colnames(all_fdr)=c('FDR_value',"Position")
row.names(all_fdr)=colnames(entropia_matrix_select)
all_fdr=all_fdr[order(all_fdr$'FDR_value'),]
return(all_fdr)
}
plot_losing_score=function(position,entropia_matrix,cluster,time,fit,indici){
position=colnames(entropia_matrix)[colnames(entropia_matrix)==position]
cluster=cluster[as.numeric(indici[[position]])]
cluster=factor(cluster,levels=c("Winner Epiblast","Intermediate","Loser Epiblast"),ordered = TRUE)
mutation=entropia_matrix[as.numeric(indici[[position]]),position]
time=time[as.numeric(indici[[position]])]
dati_scatter_1=data.frame(time,mutation)
plot<-ggplot(dati_scatter_1, aes(x=time, y=mutation)) + geom_point(aes(colour = (cluster)))+labs(col="Cluster")+labs(x="Losing Score")+labs(y="Heteroplasmy") +scale_color_manual(values=c("#DD6400", "#0000FF","#006400"),labels=c("Winning Epiblast","Intermediate", "Losing Epiblast"))+ggtitle(paste(as.character(position),"adjusted p value=",round(fit[position,1],(floor(-log10(fit[position,1])))+1),sep=" "))+
geom_smooth(method="loess",se=F,color="black",fullrange=F)
plot
}
plot_batch=function(position,entropia_matrix,batch,colour,indici,cluster,number){
position=colnames(entropia_matrix)[colnames(entropia_matrix)==position]
batch=batch[as.numeric(indici[[position]])]
mutazione=entropia_matrix[as.numeric(indici[[position]]),position]
colour=colour[as.numeric(indici[[position]])]
cluster=cluster[as.numeric(indici[[position]])]
dati_scatter_1=data.frame(batch,mutazione)
ggplot(data =dati_scatter_1
,aes(x = batch
,y = mutazione
,col =factor(cluster,levels=c("Winner Epiblast","Intermediate","Loser Epiblast"),ordered = TRUE))) +
geom_boxplot(width=0.9,alpha=0.8)+
labs(col="Cluster")+labs(x="Batch")+labs(y="Heteroplasmy") +ggtitle(position) +scale_color_manual(values=c("#DD6400", "#0000FF","#006400"),labels=c("Winner Epiblast","Intermediate", "Loser Epiblast"))+theme(axis.text=element_text(size=number))
}
plot_losing_score_average=function(average,cluster,time,indici,p_value){
cluster=cluster[as.numeric(indici)]
cluster=factor(cluster,levels=c("Winner Epiblast","Intermediate","Loser Epiblast"),ordered = TRUE)
mutation=average[as.numeric(indici)]
time=time[as.numeric(indici)]
dati_scatter_1=data.frame(time,mutation)
plot<-ggplot(dati_scatter_1, aes(x=time, y=mutation)) + geom_point(aes(colour = (cluster)))+labs(col="Cluster")+labs(x="Losing Score")+labs(y="Heteroplasmy") +scale_color_manual(values=c("#DD6400", "#0000FF","#006400"),labels=c("Winning Epiblast","Intermediate", "Losing Epiblast"))+ggtitle(paste('Average',p_value,sep=" "))+
geom_smooth(method="loess",se=F,color="black",fullrange=F)
plot
}
Only CI treated epiblast cells
entropia_matrix_fmk=entropy_method(all_pos,cells_fmk,vedere_all,vedere_numb,50,2000,filtering = 2,cell_clu_1=cell_clu_1,cell_clu_3=cell_clu_3,cell_clu_4=cell_clu_4,selection = FALSE)
entropia_matrix_ci=entropia_matrix_fmk[[1]]
allele_matrix_ci=entropia_matrix_fmk[[2]]
indici_ci=entropia_matrix_fmk[[3]]
cluster_ci=as.character(provo_4[row.names(entropia_matrix_ci),]$cluster)
cluster_ci[cluster_ci=="1"]="Winner Epiblast"
cluster_ci[cluster_ci=="3"]="Intermediate"
cluster_ci[cluster_ci=="4"]="Loser Epiblast"
condition_ci=as.character(provo_4[row.names(entropia_matrix_ci),]$condizione)
condition_ci[condition_ci=='FMK']='CI-treated'
condition_ci[condition_ci=='DMSO']='DMSO-treated'
#plot_heteroplasmy_cell_competition('327',entropia_matrix_ci,cluster_ci,indici_ci)
#plot_heteroplasmy_cell_competition('326',entropia_matrix_ci,cluster_ci,indici_ci)
#plot_allele_frequency_cell_competition('326',entropia_matrix_ci,allele_matrix_ci,cluster_ci,indici_ci#,number=5)
#plot_allele_frequency_cell_competition('327',entropia_matrix_ci,allele_matrix_ci,cluster_ci,indici_ci#,number=5)
Losing score analysis
library(gam)
#setwd("/Users/gabriele.lubatti/Desktop/Phd/Cell_Competition/Script")
#load("time_fmk_epi.Rda")
#load("cells_epi_fmk.Rda")
row.names(provo_4)=provo_4$nomi_vecchi
provo_6=provo_4[cells_epi_fmk,]
provo_6$tempo=time_fmk_epi
row.names(provo_6)=provo_6$Cellname
provo_7=provo_6[row.names(entropia_matrix_ci),]
fit_result=gam_fit_cell_competition(entropia_matrix_ci,indici_ci,provo_7$tempo,0.01,10)
fit_result_final=fit_result[fit_result$FDR_value<0.001,]
Figure 6 B, C, D, E, F, G and Extended Figure 7 A, B, C, D, E
p=list()
for(i in row.names(fit_result_final)){
q=plot_losing_score(i,entropia_matrix_ci,cluster_ci,provo_7$tempo,fit_result_final,indici_ci)
p=list(p,q)
}
p
## [[1]]
## [[1]][[1]]
## [[1]][[1]][[1]]
## [[1]][[1]][[1]][[1]]
## [[1]][[1]][[1]][[1]][[1]]
## [[1]][[1]][[1]][[1]][[1]][[1]]
## [[1]][[1]][[1]][[1]][[1]][[1]][[1]]
## [[1]][[1]][[1]][[1]][[1]][[1]][[1]][[1]]
## [[1]][[1]][[1]][[1]][[1]][[1]][[1]][[1]][[1]]
## [[1]][[1]][[1]][[1]][[1]][[1]][[1]][[1]][[1]][[1]]
## [[1]][[1]][[1]][[1]][[1]][[1]][[1]][[1]][[1]][[1]][[1]]
## list()
##
## [[1]][[1]][[1]][[1]][[1]][[1]][[1]][[1]][[1]][[1]][[2]]
## `geom_smooth()` using formula 'y ~ x'

##
##
## [[1]][[1]][[1]][[1]][[1]][[1]][[1]][[1]][[1]][[2]]
## `geom_smooth()` using formula 'y ~ x'

##
##
## [[1]][[1]][[1]][[1]][[1]][[1]][[1]][[1]][[2]]
## `geom_smooth()` using formula 'y ~ x'

##
##
## [[1]][[1]][[1]][[1]][[1]][[1]][[1]][[2]]
## `geom_smooth()` using formula 'y ~ x'

##
##
## [[1]][[1]][[1]][[1]][[1]][[1]][[2]]
## `geom_smooth()` using formula 'y ~ x'

##
##
## [[1]][[1]][[1]][[1]][[1]][[2]]
## `geom_smooth()` using formula 'y ~ x'

##
##
## [[1]][[1]][[1]][[1]][[2]]
## `geom_smooth()` using formula 'y ~ x'

##
##
## [[1]][[1]][[1]][[2]]
## `geom_smooth()` using formula 'y ~ x'

##
##
## [[1]][[1]][[2]]
## `geom_smooth()` using formula 'y ~ x'

##
##
## [[1]][[2]]
## `geom_smooth()` using formula 'y ~ x'

##
##
## [[2]]
## `geom_smooth()` using formula 'y ~ x'

Extended Figure 7 F, G, H, I, J, K
utile=provo_7
numero=rep(0,length(utile$cluster))
numero[utile$cluster==1&utile$Cell_type=="FMK 3"]="a)Win 1"
numero[utile$cluster==1&utile$Cell_type=="FMK 4"]="a)Win 5"
numero[utile$cluster==1&utile$Cell_type=="FMK 5"]="a)Win 4"
numero[utile$cluster==1&utile$Cell_type=="FMK 6"]="a)Win 2"
numero[utile$cluster==1&utile$Cell_type=="FMK 7"]="a)Win 3"
numero[utile$cluster==3&utile$Cell_type=="FMK 3"]="b)Int 1"
numero[utile$cluster==3&utile$Cell_type=="FMK 4"]="b)Int 5"
numero[utile$cluster==3&utile$Cell_type=="FMK 5"]="b)Int 4"
numero[utile$cluster==3&utile$Cell_type=="FMK 6"]="b)Int 2"
numero[utile$cluster==3&utile$Cell_type=="FMK 7"]="b)Int 3"
numero[utile$cluster==4&utile$Cell_type=="FMK 3"]="c)Los 1"
numero[utile$cluster==4&utile$Cell_type=="FMK 4"]="c)Los 5"
numero[utile$cluster==4&utile$Cell_type=="FMK 5"]="c)Los 4"
numero[utile$cluster==4&utile$Cell_type=="FMK 6"]="c)Los 2"
numero[utile$cluster==4&utile$Cell_type=="FMK 7"]="c)Los 3"
colore=rep(0,length(numero))
colore[grep("Int",numero)]="#0000FF"
colore[grep("Los",numero)]="#006400"
colore[grep("Win",numero)]="#DD6400"
plot_batch('326',entropia_matrix_ci,numero,colore,indici_ci,cluster_ci,number=6)

plot_batch('327',entropia_matrix_ci,numero,colore,indici_ci,cluster_ci,number=6)

plot_batch('303',entropia_matrix_ci,numero,colore,indici_ci,cluster_ci,number=6)

plot_batch('304',entropia_matrix_ci,numero,colore,indici_ci,cluster_ci,number=6)

plot_batch('305',entropia_matrix_ci,numero,colore,indici_ci,cluster_ci,number=6)

plot_batch('300',entropia_matrix_ci,numero,colore,indici_ci,cluster_ci,number=6)

Figure 6 I
index_rnr1=Reduce(intersect, list(as.numeric(indici_ci[[row.names(fit_result_final)[1]]]),as.numeric(indici_ci[[row.names(fit_result_final)[2]]]),as.numeric(indici_ci[[row.names(fit_result_final)[6]]]),as.numeric(indici_ci[[row.names(fit_result_final)[7]]]),as.numeric(indici_ci[[row.names(fit_result_final)[8]]]),as.numeric(indici_ci[[row.names(fit_result_final)[11]]])))
position_rnr1=c(row.names(fit_result_final)[1],row.names(fit_result_final)[2],row.names(fit_result_final)[6],row.names(fit_result_final)[7],row.names(fit_result_final)[8],row.names(fit_result_final)[11])
u=entropia_matrix_ci[index_rnr1,position_rnr1]
studio=data.frame(u)
#la parte sulla correlazione non viene per il formato del dataframe
correlazione=cor(studio,method="spearman")
colnames(correlazione)=paste((substr(colnames(correlazione),2,nchar(colnames(correlazione)))),c("mt-Rnr1","mt-Rnr1","mt-Rnr1","mt-Rnr1","mt-Rnr1","mt-Rnr1"))
row.names(correlazione)=paste((substr(row.names(correlazione),2,nchar(row.names(correlazione)))),c("mt-Rnr1","mt-Rnr1","mt-Rnr1","mt-Rnr1","mt-Rnr1","mt-Rnr1"))
ht21 = Heatmap(correlazione# heat.vals_Mutant_Paired #
,cluster_rows = TRUE
, col= colorRamp2(c(0, round(max(correlazione))), c("white", "red"))
, name = "Spearman correlation"
, column_title = "Spearman correlation FMK cells (215 cells)"
#, cluster_columns = as.dendrogram(my.tree)
, cluster_columns = TRUE
#, cluster_rows = TRUE #=as.dendrogram(as.numeric(row.names(log_data_tpm)[1:10]) )
#, top_annotation = haCol1
, row_names_gp = gpar(fontsize = 8
#,col = geneColors
)
,column_names_gp = gpar(fontsize=8)
, show_column_names = T
, show_row_names = T
)
draw(ht21
#,column_title = "Absolute values", column_title_side = "top"
)

Extended Figure 7 I
index_top=Reduce(intersect, list(as.numeric(indici_ci[[row.names(fit_result_final)[1]]]),as.numeric(indici_ci[[row.names(fit_result_final)[2]]]),as.numeric(indici_ci[[row.names(fit_result_final)[3]]]),as.numeric(indici_ci[[row.names(fit_result_final)[4]]]),as.numeric(indici_ci[[row.names(fit_result_final)[5]]]),as.numeric(indici_ci[[row.names(fit_result_final)[6]]]),as.numeric(indici_ci[[row.names(fit_result_final)[7]]]),as.numeric(indici_ci[[row.names(fit_result_final)[8]]]),as.numeric(indici_ci[[row.names(fit_result_final)[9]]]),as.numeric(indici_ci[[row.names(fit_result_final)[10]]]),as.numeric(indici_ci[[row.names(fit_result_final)[11]]])))
position_top=c(row.names(fit_result_final)[1],row.names(fit_result_final)[2],row.names(fit_result_final)[3],row.names(fit_result_final)[4],row.names(fit_result_final)[5],row.names(fit_result_final)[6],row.names(fit_result_final)[7],row.names(fit_result_final)[8],row.names(fit_result_final)[9],row.names(fit_result_final)[10],row.names(fit_result_final)[11])
u=entropia_matrix_ci[index_top,position_top]
studio=data.frame(u)
#la parte sulla correlazione non viene per il formato del dataframe
correlazione=cor(studio,method="spearman")
colnames(correlazione)=paste((substr(colnames(correlazione),2,nchar(colnames(correlazione)))),c("mt-Rnr1","mt-Rnr1","mt-Rnr1","mt-Rnr1","mt-Rnr1","mt-Rnr1"))
row.names(correlazione)=paste((substr(row.names(correlazione),2,nchar(row.names(correlazione)))),c("mt-Rnr1","mt-Rnr1","mt-Rnr1","mt-Rnr1","mt-Rnr1","mt-Rnr1"))
ht21 = Heatmap(correlazione# heat.vals_Mutant_Paired #
,cluster_rows = TRUE
, col= colorRamp2(c(0, round(max(correlazione))), c("white", "red"))
, name = "Spearman correlation"
, column_title = "Spearman correlation FMK cells (214 cells)"
#, cluster_columns = as.dendrogram(my.tree)
, cluster_columns = TRUE
#, cluster_rows = TRUE #=as.dendrogram(as.numeric(row.names(log_data_tpm)[1:10]) )
#, top_annotation = haCol1
, row_names_gp = gpar(fontsize = 8
#,col = geneColors
)
,column_names_gp = gpar(fontsize=8)
, show_column_names = T
, show_row_names = T
)
draw(ht21
#,column_title = "Absolute values", column_title_side = "top"
)
